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ABSTRACT 


This  research  is  concerned  with  two  log  spectral  estimators  in 
the  context  of  both  stationary  and  nonstationary  signals.  They 
differ  because  in  one  smoothing  is  realized  before  the  logarithmic 
transformation,  while  the  other  is  smoothed  in  the  logarithmic 
domain.  It  is  shown  that  for  stationary  signals  the  two  estimators 
are  similar,  differing  in  expected  value  by  only  a universal 
constant.  The  first  estimator,  however,  is  smoother.  For 
nonstationary  signals,  the  estimators  are  biased  by  different 
amounts  dependent  upon  the  nonstationarity.  The  difference  between 
the  estimators  is  shown  to  be  a sensitive  test  for  nonstationarity. 
The  estimators  are  used  in  the  analysis  and  implementation  of  two 
solutions  to  the  problem  of  blind  deconvolution.  It  is  found  that 
the  methods  are  equivalent  for  stationary  signals,  but  differ 
markedly  for  nonstationary  signals  in  the  presence  of  stationary 
background  noise.  Recommendations  are  made  for  the  practical 
digital  implementation  of  the  log  spectral  estimators. 
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CHAPTER  1 


INTRODUCTION  AND  OVERVIEW 

Spectral  estimation  is  a well  knoun  and  commonly  used  technique 
for  data  analysis  (e.g.,  see  til-141).  With  the  advent  of  digital 
signal  processing  (e.g.,  see  IS!  , tGl  ) ancj  the  development  of 
high-speed  techniques  such  as  the  Fast  Fourier  Transform  t71  , 
digital  algorithms  for  spectral  analysis  have  been  implemented  and 
many  new  applications  discovered.  This  research  is  concerned  with 
two  particular  estimators  of  log  spectra  and  their  application  to 
digital  signal  processing. 

rtuch  is  known  about  the  statistics  of  spectral  estimation. 
Conventional  estimators,  however,  are  often  limited  to  stationary 
(time  invariant)  signals  or,  at  most,  to  specific  types  of 
nonstationary  processes  131,  IS) . Since  many  practical  signals,  such 
39  speaking  or  singing,  exhibit  complex  nonstationarity,  these 
estimators  may  actually  be  misleading. 

Data,  including  spectral  estimates,  are  often  presented  on  a 
logarithmic  scale  since  such  representations  not  only  have  a smaller 
dynamic  range  but  frequently  a variance  independent  from  the  data. 
Log  spectral  estimates  are  often  computed  by  transforming  a spectral 
estimator,  smoothed  by  averaging,  into  the  logarithmic  domain. 
Research  has  been  published,  however,  in  which  data  (including 
spectral  estimates)  are  averaged  in  the  logarithmic  domain  (e.g.. 
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see  [9] -[14]).  As  we  shall  see,  particularly  for  nonstaticmary 
signals,  the  order  of  averaging  may  produce  significantly  different 
resul ts. 

In  this  research,  we  are  concerned  with  the  statistical 
analysis  of  two  similar  but  different  log  spectral  estimators.  The 
first,  for  convenience  termed  the  log  average  spectrum,  is  the 
logarithm  of  a conventional  smoothed  spectral  estimator.  The 
second,  the  average  log  spectrum,  differs  in  the  fact  that  smoothing 
is  done  in  the  logarithmic  domain.  Ue  consider  the  properties  of 
these  estimators  not  only  for  stationary  signals,  but  in  terms  of  a 
model  of  nonstationari ty  in  which  the  energy  at  each  frecjuency  is 
allowed  to  vary  slowly  with  time. 

A useful  application  of  these  estimators  is  in  an  area  of 
signal  processing  known  as  blind  deconvolution  115) , the  problem  of 
separating  two  convolved  signals  where  neither  is  known  a priori.  A 
knowledge  of  the  properties  of  log  spectral  estimation  is  helpful  in 
understanding  two  particular  approaches  to  this  problem. 

These  topics  are  discussed  in  two  sections.  The  first 
developes  the  statistical  analysis  of  the  estimators.  The  second 
discusses  their  digital  implementation  and  application  to  blind 
deconvolution.  In  the  first  section,  chapter  2 is  a summary  of 
fundamental  results  from  probability  and  statistics,  including  a 
discussion  of  random  processes  and  conventional  techniques  of 
spectral  estimation.  Since  this  material  is  available  from  any  of 
several  excellent  texts,  some  of  which  are  referenced  in  chapter  2, 
it  is  presented  with  a minimum  of  detail  and  mathematical  rigor.  It 
is  intended  to  establish  a common  vocabulary  and  lay  groundwork  for 
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the  remainder  of  this  work. 

Chapter  3 is  a detailed  discussion  of  the  two  log  spectral 
estimators  in  terms  of  stationary,  Gaussian  processes.  Included  are 
derivations  of  pertinent  statistics.  Empirical  results  for  computer 
generated  signals  are  presented. 

Chapter  4 contains  an  extension  of  these  results  into  the 
domain  of  nonstationarity.  A simple  model  is  proposed  and  the 
statistical  properties  of  the  estimators  derived.  Empirical  results 
are  presented  for  some  practical  signals  as  well  as  for  simulated 
data. 

In  the  second  section,  chapter  5 presents  an  application  of 
these  results  to  digital  signal  processing  in  the  context  of  blind 
deconvolution.  The  problem  is  discussed  and  two  solutions  are 
outlined:  one  based  on  the  homomorphic  filtering  theory  of 
Oppenheim,  et  al.  1161,  and  the  other  on  the  application  of 
conventional  spectral  estimates.  An  analysis  of  these  two  solutions 
is  discussed  in  terms  cf  log  spectral  estimates.  The  impact  of 
nonstationari ty  and  additive,  stationary  noise  is  included,  and  the 
results  of  deresonating  an  old  acoustic  recording  of  the  famous 
tenor  Enrico  Caruso  are  compared  to  an  experiment  simulating  the 
actual  data. 

Finally,  chapter  G presents  a summary  of  the  conclusions  of 
this  research.  A part  of  this  summary  includes  suggested  procedures 
for  practical  digital  spectral  estimation. 


I 


CHAPTER  2 

STATISTICAL  FUNDAMENTALS 

2.1  Random  Variables  anc'  Probability  Distributions 

For  a given  experiment,  the  set  of  all  possible  outcomes  is 
called  a sample  space.  These  outcomes  may  be  grouped  in  various 
nays  to  form  events  which  have  a probability  of  occurance  between 
zero  and  one  [3,pp.5S-57l  . A random  variable,  X(k),  is  a function 
associating  a real  number  between  -»  and  +w  with  each  outcome,  k,  in 
the  sample  space.  In  general,  a random  variable  may  be  defined  for 
either  a continuous  or  discrete  sample  space. 

For  a discrete  random  variable,  we  define  the  probability 
distribution  fx(x)  as  the  probability  that  the  random  variable  X 
takes  on  the  particular  value  x.  The  distribution  fx(x)  has  the 
properties  C17,p.l55] 

fx(x)  > 0 for  ail  x (2.1a) 

and 

£fx(x)  = 1 (2.1b) 

where  the  summation  in  (2.1b)  is  over  all  x. 

For  a continuous  random  variable,  it  becomes  meaningless  to 
talk  about  such  a frequency  distribution  function.  However,  we  can 
define  a cumulative  distribution  function  (cdf),  Fx(x),  where  Fx(x) 
represents  the  probability  that  the  random  variable  X has  a value 
less  than  or  equal  to  x l4,p.B2),  i.e.. 
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Fx(x)  ■ Prob  IX  < x)  . (2.2) 

Note  that  as  x -•  Fx(x)  -•  0 and  as  x -•  +»,  F„(x)  -*  1. 

If  the  cdf  is  smooth  enough  to  be  differentiable,  we  can  also 
define  a probability  density  function  (pdf)  f„(x)  for  a continuous 
random  variable  (4,p.B2) 

fx(x)  = dFx(x)/dx,  for  almost  all  x.  (2.3) 

Althougn  not  a distribution,  the  pdf  of  X may  be  used  to  calculate 
the  probability  that  x,  < X < x,  by  noting 

Prob  (x,  < X < x?)  » J'"’fx(x)dx  . (2.4) 

A pdf  has  properties  similar  to  those  of  the  probability 


distribution  for  discrete  random  variables  (4,p.63),  namely 

fx(x)  > 0 for  all  x (2.5a) 

and 

X:fx(x)dx  = 1 . , (2.5b) 

1 f V is  a function  of  the  random  variable  X,  i.e.,  V = g(X), 
and  g(X)  is  one-to-one,  differentiable,  and  either  monotonical  ly 
increasing  or  decreasing  then  fH(y)  is  related  to  fx(x) 
by  [17, p.312J 

fw(y)  - f* (g‘l  (y) ) -Id g‘l  (y)  / dyl  . (2.B) 


There  are  several  probability  density  functions  that  are  common 
and  particularly  useful  to  this  research.  One  of  the  most  important 
is  the  normal  or  Gaussian  distribution  denoted  N ( // . ct-1  ) and  given 
by  Q7,  p.220) 

f x ( x ) - (2nerJ) ",/J.  exp  t- (x  - lx|  < «.  (2.7) 

See  figure  2.1(a).  The  normal  distribution  is  completely  specified 
by  the  parameters  A and  r.  I t is  particularly  important  because  of 
the  Central  Limit  Theorem  which  states,  in  one  form,  that  sums  of 


independent,  identically  distributed  random  variables  with  finite 
means  and  variances  quickly  tend  to  a normal  distribution  regardless 
of  the  initial  distribution  [17, p. 4311.  Consequently,  it  is  often 
possible  to  describe  simple  sums  (such  as  averages)  with  Gaussian 
statistics. 

Some  other  useful  distributions  are  the  uniform  [17, p. 2201, 
f x (x)  = 1 / ( b - a),  a < x < b 

= 0,  otherwise,  (2.8) 

the  lognormal  [18, p. 8), 

f x ( x ) = (xV*2n)  ‘1/J>exp  C-  ( log  x - /u ) z / 2cr2]  , x > 0 

= 0,  otherwise  (2.9) 

(If  X = exp  (Y)  and  V ~ N ( // , crz ) . then  X has  a normal  distribution; 
see  figure  2.1(b).)  and  the  exponential  [17, p. 220) 
f x ( x ) ■*  (1  / fj) -exp  (-x  / /d  , x > 0 

= 0 , otherwise  (2.10) 

(see  figure  2. 1(c)). 

Let  Z,,  Z ,,  •••,  Z„  be  mutually  independent  normal  random 

variables  such  that  Z,  ~ N(0,1).  Then  X = Z,2  + Z,2  + •••  + Zn2 

has  a chi-square  distribution  [4,p.79]  with  n degrees  of  freedom 

given  by 

fy(x)  = [2"'2r(n/2)]-,(x)",2',-ev.p(-x/2) , x > 0 
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0,  otherwise 

(2.11) 

where  T(x)  is 

the 

gamma 

function  (see  Appendix  A). 

Substitution  of 

n » 2 shows 

that 

the 

exponential  distribution 

is  a 

chi-square 

distribution 

wi  th 

V*  m 

2 (see  figure  2.1(c)). 

The 

chi-square 

distribution 

is  i 

t se  1 f 

a special  case  of  the  more  general  gamma 

distribution 

117.  p 

.181) . 

If  X > Z,2  + Z?2  + . — 

+ Z 

„2  and  Z,  ~ 

7 

NdD.o-*)  with  er*  not  necessarily  1,  then  X/o-  is  a chi-square  random 
variable  with  n degrees  of  freedom  (denoted  x?n) . If  the  Z' s are 
not  zero  mean,  then  a new  distribution  arises  termed  a non-central 
chi-square  distribution  119, p. 5441. 

If  X - r<X2„,  then  Y =*  logX  is  distributed  as  log  chi-square 
with  n degrees  of  freedom  [20, p. 251  (Y  * log(r*x2n)).  As  shown  in 
Appendix  B,  the  log  chi-square  pdf  has  the  form 
fw(y)  = T"'  (n  / 2)  [exp  (y  - log2r] 

exply  - log2r  - exp(y  - log2r)]  . (2.12) 

For  n - 2,  (2.12)  reduces  to  the  interesting  form  (figure  2.1(d)) 

f a ( y ) = exp  [y  - logyi/x  - exp(y  - log//*))  (2.13) 

where  ■>  2r  is  the  expectation  of  X (see  section  2.2). 

When  more  than  one  random  variable  is  defined  on  a sample 
space,  a joint  probabi  1 i ty  density  function  may  be  defined.  The 
joint  cdf  and  pdf  of  two  random  variables,  X and  Y,  are  denoted 
Fxv^x.y)  and  fx«(x,y),  respectively,  and  given  by  [4,p.B51 

Fxy(x,y)  = Pro))  IX  < x,  Y < y)  (2.14a) 

and 

fx„(x,  y)  = 3JFXJ(x,  y)  / 3x3y.  (2.14b) 

If  f**(x,y)  - fx(x)>fy(y)  then  X and  Y are  said  to  be  statistically 
independent  [17, p . 295] . 

2.2  Statistical  Parameters 

Although  a random  variable  is  described  by  its  probability 
density  function  (or  cdf),  it  is  useful  to  define  various  parameters 
which  help  characterize  it.  Among  the  most  common  of  these  are  the 
mean  and  variance.  Suet  simple  parameters  are  particularly  useful 
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in  charac terizing  data  when  the  pdf  is  not  explicitly  Known. 
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FIGURE  2.1 


Selected  probability  density  functions:  (a)  Gaussian,  (b)  log- 
normal, (c)  exponential  or  chi-square  with  2 degrees  of  freedom,  and 
(d)  log  chi-square  with  2 degrees  of  freedom. 
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The  expected  value  lor  expectation)  of  a function,  g(X),  of  the 
continuous  random  variable  X is  defined  as 

E(g(X))  - X“g  (x)  fx  (x)  dx  . (2.15) 

For  a discrete  random  variable,  (2.15)  becomes 

E (g  (X) ) - 2g(x)fx(x)  (2.1G) 

where  the  summation  in  (2.16)  is  over  all  x.  The  expectation 
operator  has  the  important  property  of  being  linear.  I f X,  V are 
random  variables  and  r,,  r2  are  constants,  then 

E(r,X  + r7Y!  - r,E  (X)  + r7E(Y)  (2.17) 

for  all  X,  Y,  rt,  and  r?  '.17, p. 355). 

If  g (X)  = Xn,  then  Elg(X))  = E SX")  = X:xnfx(x)dx  is  called  the 
nts  moment  of  X.  For  n = 1 we  have  n = E (X)  = X'^fxlxldx  where  n is 
the  average  or  mean  value  of  the  random  variable  X.+  For  n = 2, 
E (X2)  is  the  mean  square  value. 

Proceeding  similarly,  E((X  - Ax)n)  is  the  nth  central  moment  of 
X.  For  n - 2,  we  define  the  variance  (o-2)  as  var  IX)  * E f (X  - //)2) 
and  the  standard  deviation  (o-)  as  the  positive  square  root  of  the 
variance.  It  is  easily  shown  that  var  IX)  = E IX2)  - E2  !X) . 

Qualitatively,  the  mean  is  a simple  average,  the  mean  square  a 
measure  of  the  general  intensity  of  the  random  variable,  and  the 
standard  deviation  a measure  of  the  spread  about  the  mean.  Other 
parameters  often  given  fcr  a random  variable  include  the  median  (the 
value  of  x for  which  there  is  equal  area  on  either  side  under  the 
pdf),  and  the  mode  (a  value  of  x for  which  the  pdf  has  a relative 


tThis  and  following  results  are  given  for  continuous  random 
variables  only.  Similar  results  apply  to  discrete  random  variables 
with  appropriate  summations  replacing  the  integrals. 
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maximum)  [17, pp. 205, 213)  . 

It  is  frequently  useful  to  characterize  a random  variable,  X, 
as  being  proportional  or  approximately  proportional  to  a chi-square 
random  variable  with  n degrees  of  freedom,  i.e.,  X = r*xz„.  This  is 
particularly  useful  in  computing  the  statistical  properties  of 
spectral  estimates  [21 , pp. 21-251  , [221 . To  do  so,  a useful  concept. 
Sat  ter thwai te’ s approximation,  is  used  to  estimate  r and  n. 
Specifically 

n - EOF  (X)  - 2-EMX)  /var  (X)  . (2.18a) 

and 

r - var  (X)  / (2-EIX) ) =EfXl/n  (2.18b) 

where  EDF  )X)  is  the  equivalent  degrees  of  freedom  of  X [8,  p. 2731. 
This  approximation  is  used  in  chapter  4 when  developing  the 
statistics  of  log  spectral  estimators  for  nonstationary  processes. 

The  moments  of  a given  pdf  often  have  a simple  functional 
relationship  to  the  pa-ameters  characterizing  the  distribution. 
Table  2.1  lists  the  mean  and  variance  for  the  pdfs  given  in  section 
2.1  with  reference  to  the  indicated  equations  [17, pp. 220, 3483  . 

These  definitions  can  be  extended  to  two  or  more  random 
variables.  Thus,  as  an  example,  for  two  random  variables  X and  Y, 
we  define  the  product  moment  and  covariance  as 

EIX.Y)  - X:X:x'yf«(x.y)dxdy  (2.19a) 

and 

coviX.YI  - E ( (X  - • <Y  - a*>)  (2.19b) 

where  fXM(x,y)  is  the  joint  pdf  of  X and  Y and  Ax  and  Ax  sre  their 


respective  means.  It  can  be  shown  [17, p. 3553  that 
var  (X  + Y)  - var(X)  + variY)  + 2-covtX.Y) 


(2.20a) 
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and 

cov  (X,  Y)  = E<X-Y)  - E (X)  E (Yf  . (2.20b) 

If  covlX.Yl  = 0,  then  X and  Y are  said  to  be  linearly 
independent  (or  uncorrelated).  If  X and  Y are  statistically 
independent,  then  E (X*Y)  = EfXlElYl  (since  fxU(x,y)  = fx!*)fM(y)) 
and,  from  (2.20b),  cov(X,YI  = 0.  Thus,  X and  Y are  also 
uncorrelated.  Note,  however,  that  except  under  certain  conditions 
(e.g.,  normality),  the  reverse  is  not  necessarily  true  C4,p.74). 


TABLE  2.1 

MEAN  ANO  VARIANCE  FOR  SELECTEO 
PROBABILITY  DENSITY  FUNCTIONS 


PDF 

Clean 

Variance 

Normal  (2.7) 

A 

<72 

Uniform  (2.8) 

(b  + a)  / 2 

(b  - a) 2 / 12 

Lognormal  (2.9) 

exp(/i  + l/2cr2) 

exp  {2m  + 2 cr2)  - 

exp  (2 m + <7'J) 

Exponential  (2.10) 

4 

A1 

r.x’n  (2.11) 

r-n 

2r2n 

log  (r*x7„)  (2.12)  + 

V'  (n  / 2)  + 

log (2r) 

V''  (n  / 2) 

log  (r*x??)  (2.13)  + 

log2r  - y 

n2/G 

2.3  Statistical  Estimates 

In  experimental  work,  it 

is  often  desirable 

to  estimate  from  a 

+ V - Euler’s  constant  (v  - 

0.57721  — ) and 

\M  t) 

, (A'lt)  are  the 

digamma  and  trigamma  functions 

, respectively 

( see 

Appendix  A) . 
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given  set  of  data  (representing  individual  sample  values  of  a random 
variable)  one  or  more  of  the  associated  statistical  parameters.  The 
development  and  analysis  of  suitable  estimators  is  the  domain  of 
statistical  estimation  theory. 

Since  a statistical  estimate  is  derived  from  random  samples,  it 
is  itself  a random  variable  with  an  associated  pdf  (called  the 
sampling  distribution),  riean  value,  variance,  etc.  If  a'  estimates 
a,  then  a'  is  said  to  be  consistent  if  E((a'  - a)*)  -•  0 as  the 
number  of  samples  increases  and  unbiased  if  the  expected  value  of  ex' 
• equals  a [3,  pp.  100-11  . For  example,  it  can  be  shown  that  m *= 

(l/N)  2iNi*i«  the  sample  mean,  is  both  consistent  and  unbiased  since 
E(m)  - fx  and  El(m  - a)2)  - 0 as  N " N is  the  number  of  samples, 

Xt- 

I 

In  choosing  an  estimator  for  a particular  parameter  ot,  it  is 
desirable  to  minimize  the  variance  and  the  bias 

B fa ' ) = Efa'l  - a (2.21) 

where  8(a')is  the  bias,  a the  parameter  and  a'  an  estimator  of  a. 
Unfortunately,  it  is  often  the  case  that  an  unbiased  estimator  is 
not  the  one  with  the  smallest  variance,  and  vice  versa. 

Accordingly,  several  different  criteria  have  been  devised  for 

describing  an  estimator  and  compromising  between  these  two 

conditions.  For  example,  the  mean  square  error  [4, p.983  can  be 
minimized  where  msefa')  « El(a'  - a)2)  = varfa'l  + BMa'J.  Another 
method  is  to  compute  the  likelihood  function  [4, pp. 99-1021 , and  use 
the  maximum  of  this  function  as  the  desired  estimate.  In  this  work, 
we  discuss  both  the  bias  and  variance  of  the  log  spectral  estimators 
and  are  not  generally  concerned  with  maximum  likelihood 
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considerations. 

In  the  experimental  aspects  of  this  research,  the  mean  and 
variance  of  experimental  data  are  frequently  estimated.  The 
formulas  used  for  these  computations  are  standard  and  given 
by  [3, pp. 100-21 

m - (l/N)ZiNx><i  (2.22a) 

and 

s1  - [1/  (N-niZi"  <X;  - m) 2 (2.22b) 

where  m is  the  sample  mean,  sJ  the  sample  variance  and  N the  number 
of  samples,  xt.  The  factor  1/  (N— 1 ) in  (2.22b)  insures  that  the 
sample  variance  is  unbiased  (although  it  does  not  have  the  smallest 
variance) . 

The  above  discussions  concern  point  estimates,  i.e.,  a single 
value  as  an  estimate  for  a statistical  parameter.  It  is  frequently 
more  useful  to  generate  a confidence  interval  and  specify  the 
probability  that  the  desired  parameter  falls  within  that  interval. 
Computation  of  confidence  intervals  requires  that  the  pdf  of  the 
sample  data  be  known,  which  is  often  not  the  case.  However,  as  many 
estimates  involve  sums  of  data  (such  as  (2.22)),  application  of  the 
Central  Limit  Theorem  enables  us  to  assume  a normal  distribution  for 
the  point  estimates  and  derive  approximate  intervals.  In  several 
instances  in  this  research,  confidence  intervals  are  computed  along 
with  point  estimates.  In  doing  so,  the  following  formulas  are 
utilized  123,  pp. 296-981  . Note  that  N is  the  number  of  samples  and 
n - N - 1. 

The  100(1  - <u)%  confidence  interval  for  a sample  mean,  m,  with 
unknown  variance  is  giver  by 


tm  - t (cr  / 2;  n) -s  / Nv2  < u < m + t (a/ 2;  n)  *s  / N1/2) 


(2.23) 


where  t(a/2;n)  is  the  100(a/2)  percentage  point  for  Student’s 
t-distribution  (17, p. 180)  with  n degrees  of  freedom.  Similarly,  the 
100(1  - a)%  interval  for  a sample  variance,  s2,  is  given  by 

(n*  s2 / X?(a  / 2;  n)  < <r2  < n.S2/x2(l  - ct/2;n)]  (2.24) 

where  x7(a/2;n)  is  the  100(a/2)  percentage  point  for  the  chi-square 
distribution  with  n degrees  of  freedom.  See  (23)  for  tables  of 
percentage  points. 

Closely  allied  to  confidence  intervals  is  the  concept  of 
hypothesis  testing.  Basically,  the  idea  is  to  formulate  a 
hypothesis  regarding  a set  of  random  data,  compute  a test  statistic 
(as  a function  of  the  data,  hypothesis,  sample  size,  etc.)  and 
determine  a region  of  acceptance.  If  the  test  statistic  falls 
within  the  region,  the  hypothesis  is  accepted;  otherwise  it  is 
rejected.  In  later  discussions,  we  will  use  an  hypothesis  test 

called  a chi-square  goneness  of  fit  test  to  determine  if  observed 
data  obeys  a predicted  probability  distribution  (23,  pp.  458-GO) . 
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2.4  Random  Processes. 

Physical  phenomena  are  frequently  represented  by  a series  of 
observed  data.  This  data,  whether  continuous  (analog)  or  discrete 
(digital),  may  be  characterized  by  a functional  relationship  between 
the  independent  variable (s)  + of  the  observation  and  data  values. 


+For  a one-dimensional  process,  the  independent  variable  is  often 
time.  Except  where  ncted,  these  discussions  are  for  functions  of 
time  only.  However,  many  of  these  results  are  applicable  in  two  or 
more  dimensions  (e.g.,  see  (24)  and  (25)). 
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Data  resulting  from  one  particular  set  of  measurements  is  called  a 
sample  function;  the  collection  of  all  possible  sample  functions 
forms  a process  or  time  series. 

A time  series  is  broadly  classified  as  deterministic  or  random 
(non-deterministic) . A deterministic  process  is  one  in  which  there 
is  an  explicit  (although  not  necessarily  simple)  mathematical 
relationship  between  the  data  and  the  independent  variable.  All 
sample  functions  for  such  a process  are  determined  by  the  same 
functional  relationship.  Within  the  limits  of  the  measuring  device, 
measurements  taken  at  one  time  are  equivalent  to  those  taken  later. 
Random  processes  (or  stochastic  processes)  do  not  have  such  explicit 
functional  relationships,  but  are  sets  of  random  variables  and  are 
best  described  statistics  1 ly  [3,pp.l-14]. 

Many  common  phenomena  are  deterministic.  A vibrating  string, 
the  orbit  of  a satellite,  and  the  current  flowing  in  an  electronic 
circuit,  for  example,  are  conveniently  described  by  an  explicit 
mathematical  function  (such  as  a sinusoid  for  a simply  vibrating 
string).  Conversely,  the  location  of  an  electron  in  an  atom, 
thermal  noise  in  an  amplifier  or  the  amount  of  water  passing  a 
particular  point  in  a mountain  stream  are  examples  of  random 
processes. 

Often  the  categorization  of  physical  data  is  arbitrary  and  a 
matter  of  convenience.  One  might  argue,  for  example,  that  the 
motion  of  a vibrating  spring  is  not  exactly  sinusoidal  because  of 
random  interactions  with  molecules  of  air.  Similary,  given  enough 
information  about  the  piysical  nature  of  the  stream  bed  and  its 
surrounding  environment,  etc.,  the  flow  of  water  might  be  accurately 
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predicted  for  all  time.  Clearly,  however,  for  most  applications,  a 
simple  vibrating  string  is  best  thought  of  as  deterministic  and  the 
stream  as  random. 

These  considerations  motivate  our  classifying  signals  such  as 
accoustic  waves,  light  being  recorded  on  film  and  seismic  waves  as 
random  processes.  A particular  recording,  photograph  or  seismic 
chart  is  a finite  i~eal ization  or  sample  function  of  the  process. 

Mathematically,  a random  process  is  an  ordered  or  paramaterized 
set  of  random  variables  and  is  denoted  1 x j ( t ) J , -»  < t < »,  where  j 
indexes  sample  functions  of  the  process  [4,pp.l441.  A particular 
realization  of  the  process  is  x(t).  Frequently,  we  drop  the  braces 
and  subscript  and  simply  write  x(t)  for  both  the  process  and  a 
realization  allowing  context  to  discriminate  between  them.  For  each 
value  t,  we  associate  a random  variable  with  the  process. 

As  mentioned,  a process  may  be  either  discrete  or  continuous  in 
either  the  independent  or  dependent  variables.  Throughout  these 
discussions,  we  generally  deal  with  continuous  processes.  However, 

for — practical  application  in  digital  signal  processing,  the 

processes  must  be  discrete  in  both  the  dependent  (quantized)  and 
independent  (sampled)  variables.  With  proper  consideration  given  to 
the  problems  introduced  by  sampling  and  quantizing,  most  of  the 
analysis  may  be  directly  applied  to  discrete  as  well  as  continuous 
processes.  A discrete  process  is  denoted  (xj(n)l.  However,  to 
again  simplify  notation,  the  subscript  j is  dropped  and  the  process 
indexed  as  x (n) . 

The  collection  of  all  sample  functions  that  might  be  produced 
by  the  same  random  phenomenon  is  called  an  ensemble  13, p. 101.  While 
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each  member  of  the  ensemble  is  an  explicitly  different  time  series, 
they  all  have  the  same  statistical  properties. 

Except  for  ergodic  processes  (section  2.G)  a particular  sample 
function,  x(t),  does  not  suitably  represent  the  entire  process 
(xj(t)l.  For  a given  ensemble,  there  are  (usually)  an  infinite 
number  of  possible  sample  functions.  However,  once  a particular 
sample  function  is  realized,  it  becomes  a deterministic  function 
over  the  domain  for  which  it  is  realized,  and  may  be  treated  as 
such.  For  a particular  value  of  t,  then,  it  represents  one  sample 
value  of  the  random  variable  associated  with  that  value  of  t.  Uith 
other  sample  functions  it  may  be  used  to  estimate  the  statistics  of 
the  process  at  that  point. 

A random  process  may  be  described  to  a first  order  by  the 
probability  density  function,  f>(x,t),  associated  with  each 
paramaterized  random  variable,  x(t).  Consequently,  we  may  define 
the  simple  moments  of  the  process  at  each  time,  t,  e.g.,  the  mean, 
variance,  etc.  A process  is  more  completely  described  in  terms  of 
higher  order  statistics  by  defining  the  joint  pdf  associating  the 
random  variables  at  arbitrary  times  t,,  t,,  •••,  tn  and  the 
corresponding  multivariate  moments  [4,p.l46).  However,  as  these 
pdfs  and  moments  may  be  rather  complicated,  it  usually  suffices  to 
describe  the  simpler  first  and  second  order  moments  such  as  the  mean 
and  covariance  (section  2.7). 

A random  process  is  often  refered  to  in  terms  of  the 
statistical  properties  that  characterize  it.  For  example,  if  the 


random  variables  associated  with  the  process  obey  a normal 
distribution,  the  process  is  called  a normal  or  Gaussian  random 
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process.  Similarly,  if  the  spectrum  (defined  in  section  2.8)  is 
flat,  it  is  refered  to  as  white.  Such  descriptive  terms  will  be 
useful  in  later  discussions. 

2.5  Ensemble  and  Time  Averages 

As  mentioned  in  the  previous  section,  since  a random  process  is 
a parameterized  set  of  random  variables,  the  simple  moments  of  a 
process  may  be  computed  by  statistical  averages  called  ensemble 
averages  (5, p.382).  Accordingly,  the  expected  value  of  a process  is 
Elx(t))  = J!"xfx(x,  t)  dx  (2.25) 

where  f*(x,t)  is  the  first  order  pdf  associated  with  the  random 
variable  at  time  t.  Since  f„(x,t)  is  a function  of  time,  so  is  the 
expected  value.  Other  ersemble  averages  may  be  defined  by  extending 
(2.25) 

Eig(x(t)))  = x'“g  (x)  f>  (x,  t)  dx  . (2.26) 

In  working  with  random  processes,  it  is  often  the  case  that 
only  one  member  of  the  ersemble  is  available  making  it  impossible  to 
compute  ensemble  averages.  For  this  reason  random  processes  are 
also  characterized  by  averages  computed  over  time.  For  example,  the 
time  average  representing  the  expected  value  is 

n-  = <x(t)>  = Uff  (1  / 2T)  X;  x ( t ) dt  (2.27) 

and  the  variance  is 

crIj  = < ( x ( t ) - ny>  = liff  (1  / 2T)  S.Tr  <x  ( t ) - fJj)  * d t (2.28) 
where  the  subscript  j on  and  cr2;  indicates  that  the  averages  are 
now  a function  of  the  j1h  member  of  the  ensemble  rather  than  time. 
Another  important  time  average,  the  mean  square  value  <x*  (t)>,  is 
interpreted  as  the  average  power  in  the  process.  Other  time 
averages  are  similarly  defined  15, p. 386). 
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2.G  Stationarity  and  Ercodicity. 

In  general,  for  a given  process,  time  averages  are  not  equal  to 
ensemble  averages.  Processes  for  which  they  are  equa 1 are  said  to 
be  ergodic.  If  this  is  true  for  all  possible  statistics,  then  the 
process  is  strongly  or  strict  sense  ergodic;  if  true  for  only 
selected  statistics  then  it  is  weakly  ergodic.  For  such  processes, 
estimates  of  time  averages,  which  can  be  computed  from  one  sample 
process,  serve  well  as  estimates  of  the  corresponding  ensemble 
averages.  Ue  have  found  it  particularly  useful  in  this  research  to 
frequently  assume  a process  is  ergodic  to  enable  such  estimates  to 
be  meaningful. 

Random  processes  are  also  classified  as  stationary  or 

nonstationary.  Broadly,  a stationary  process  is  one  for  which  the 
statistical  properties  are  independent  of  time,  e.g.,  fx(x, t,)  = 
fx(x,t?)  for  all  t,,  t7.  If  all  possible  statistics  are  independent 

of  time,  then  the  process  is  strongly  or  strict  sense  stationary. 
However,  if  only  the  first  k moments  are  time  independent,  then  it 
is  weakly  stationary  to  the  k1h  order.  Note  that  while  an  ergodic 
process  is  necessarily  stationary,  the  reverse  may  not  be 

true  [3,p.891.  Interestingly,  if  a normal  process  is  weakly 
stationary,  then  it  is  also  strongly  stationary  and 
ergodic  [4, p. 1491.  In  these  discussions,  by  a stationary  process  we 
mean  one  which  is  strongly  stationary. 

Some  processes  exhibit  simple  nonstationarity.  For  example, 
the  mean  value  may  increase  by  a simple  linear  trend.  Such 

nonstationarities  may  be  easily  recognized  and  removed.  Many 

signals,  however,  such  as  the  acoustic  waves  pertinent  to  this 
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research,  exhibit  much  more  complex  nonstationarity.  Not  only  do 
the  mean  and  variance  change  with  time,  but  the  energies  associated 
with  a particular  frequency  change  from  moment  to  moment  (if  in  fact 


such  a concept  is 

preserved) . 

The  effects  of 

this 

type  of 

nonstationarity  are 

a focal  point 

of  this  research. 

As 

will  be 

seen,  such  processes  are  often  modeled  by  assuming  stationarity  over 
a short  interval. 

2.7  Autocorrelation  and  Autocovariance  Functions 

Two  important  statistical  parameters  of  a random  process  are 
the  autocorrelation  and  autocovariance  functions.  For  the  process 
x(t),  the  autocorrelation  function  is  defined  as 

Rxx(r,t)  - E{x(t)‘x(t  + r)  i (2.291 

and  the  autocovariance  function  as 

Cxx(r,t)  = E ( C ><  ( t ) - ^(t)]-[x(t  + r)  - /^(t  + r)  ] 1 (2.30) 

where  // ( t ) is  the  expected  value  of  x(t)  at  time  t and  r is  the 

displacement  or  lag.  For  zero  mean  processes,  the  covariance  and 
correlation  functions  are  equal. + For  two  processes,  x(t)  and  y(t), 
a cross-correlation  function  is  defined  as 

Rxy ( r,  1 1 = E(x;t)*y(t  + r)  1 . (2.31) 

If  x(t)  is  stationary,  then  Rxx(r,  t)  and  Cxx(r,  t)  are  functions 
of  r only,  i.e.,  Ryx(r,  t)  = Rxx(t)  and  Cxx ( r.  t ) = Cxx(r).  They  are 

both  even  functions,  since  Rxx(-r)  = Rxx(r),  Cxx(-r)  = Cxx(r),  and 

it  can  be  shown  that 

varlx(t))  = Cx<(0),  (2.32a) 


■♦"Throughout  this  research,  we  make  the  simplifying  assumption  of 
zero  mean  processes  and,  thus,  Rxx(t,  t)  = Cxx(r,  t). 
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|Rxx(t)|  s |Rkx UD)  I , (2.32b) 

and 

ICxx(rr)  I < ICXX((D)I  . (2.32c) 

A normalized  form  of  the  autocorrelation  function,  the  correlation 
coefficient  is  defined  as  p(r)  = Rxx(r) /Rxx(0) . It  has  the  property 
of  ranging  between  -1  anc  +1  [3, pp. 70-71] . 

These  functions  all  have  the  property  of  providing  a measure  of 
the  linear  dependence  between  two  processes  (cross-correlation)  or 
of  a process  with  itself  (autocorrelation)  for  a given  lag,  r.  For 
a completely  uncorrelated  process  (often  called  a purely  random 
process  or  white  noise)  Cxx(r)  = a2'§(r)  where  S(t)  is  the  Dirac 
delta  function  [4,p.l57].+ 

Two  common  estimators  of  the  autocovariance  function  (often 
called  sample  autocovarience  functions)  are 
cxx  (r)  - ( 1 / T ) Jo IT'  C x ( t ) - //„]• 

[x(t  + Irl)  - pty]  dt  (2.33a) 

and 

g'v*(t)- ■■  (1  / (T  *- -Jill)  - M»]  ’ 

(x(t  + Irl)  - piy]  dt  (2.33b) 

where  x(t)  is  a sample  function  of  the  ergodic  process  (x(t)).  It 
can  be  shown  [4,p.l75]  that 

E (cxx (t)  ) = (1  - ( Irl  / T)  ] Cxx (t)  (2.34a) 

and 

E(c'xx(r)l  - C>x (ri  . (2.34b) 

+ The  definition  of  Ihe  autocovariance  and  autocorrelation 
functions  and  the  correlation  coefficient  may  differ  somewhat  from 
author  to  author  (e.g.,  see  [3,pp.68-9],  [4, pp. 154-157] ) . 


A 
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Hence  c'xxir)  is  an  unbiased  estimator  of  Cxx(r)  (although  it  does 
not  have  the  smallest  mean  square  error)  whereas  c,x(t)  is 
asymptotically  unbiased.  The  variances  of  (2.33)  are  somewhat  more 
involved  than  the  expectations,  however  they  are  generally 

proportional  to  (1/T);  accordingly,  these  estimators  are  consistent 
since  tvar  Ic^tr)  > ] = U&’  (E  f (c'xx(r)  - cxx(t))j}]  = 0 [4 , pp.  175-8]  . 

2.8  Spectral  Density  Function 

The  Fourier  transform  of  a deterministic  function  gives  a 
frequency  distribution  of  signal  strength  (4,p.25].  However,  since 
a sample  function  of  a random  process  generally  has  infinite  energy, 
i.e..  X:x2(t)dt  = »,  its  Fourier  Transform  may  not  exist  (2G,p.4S5]  . 

The  mean  square  value  or  average  power  <xMt)>,  however,  is  finite 
(since  trx  is  finite)  and  has  a frequency  distribution  called  the 
power  spectral  density  function  or  simply  the  spectrum. 

The  spectrum,  Gx(f),  of  a stationary  random  process.  x(t),  is 
actually  defined  as  the  one-dimensional  Fourier  transform  of  the 
autocovariance  function,  t Cx,(t),  associated  with  x(t) 

Gx(  f ) * f{Cxx(-.-)l 

« X-Cyx(r)  -exp  (-2n  j f r)  dr,  Ifl  < » (2.35) 

where  f represents  frequency  and  j « ( -1 ) ,/2.  Gx(f)  shows  how  the 
variance  or  average  power  of  a process  is  distributed  with 
frequency  [4,p.217]. 

By  rewriting  (2.35)  in  terms  of  an  inverse  Fourier  transform, 

i 

+A1 ternatively , Gx(f)  is  frequently  defined  in  terms  of  the 
autocorrelation  function,  Ryx(r)  [3,p.761.  However,  this  may 
introduce  impulses  into  the  spectrum  if  E'<x(t)l  * 0.  For  zero  mean 
processes  the  two  approaches  are  equivalent. 
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we  have 

Cxx(r)  - X»Gx ( f ) - exp  (2n  j f r)  df  (2.3Ba) 

and 

Cxx (0)  = cr2x  = XrGx(f)df  . (2.36b) 

For  a purely  random  process,  Cxx(r)  = <7-VS(t)  and,  thus,  Gx(f)  = 
erJx  is  constant;  hence  the  descriptive  term  white  noise. 

The  spectrum  of  a real  valued  process  is  non-negative  and  even, 

i.  e. , 

Gx(f)  > 0 for  all  f (2.37a) 

and 

Gx(-f)  - Gx(f)  for  all  f.  (2.37b) 

Another  important  property  is  the  relationship  of  Gx(f)  to  linear 
stationary  systems.  If  H(f)  is  the  frequency  response  of  a linear 
system,  x(t)  the  input  and  y(t)  the  output,  then 

GM(  f ) = |H(f)  l2-Gx(f)  . (2.38) 

For  two  processes,  *(t)  and  y(t),  we  define  the  cross-spectrum, 
Gxv  ( f ) , as 

Gxa  ( f ) - Tf  (Cya(r) ) (2.39) 

where  Cx„(r)  is  the  cross-covariance  of  x ( t ) and  ylt)  (5,  pp.  390-94]  . 
Similar  properties  to  (2.37)  can  be  derived  for  the  cross-spec trum. 

The  preceeding  definitions  are  for  stationary  processes.  If 
x ( t ) is  nonstationary,  then  the  autocovariance  is  a function  of  two 
variables  and  a simple  spectrum  is  not  defined.  There  are  two  basic 
approaches  to  defining  tte  spectrum  of  a nonstationary  process.  One 
is  to  do  it  in  terms  of  a two-dimensional  Fourier  transform.  Being 
a function  of  two  frequencies,  however,  such  a spectrum  may  be 
difficult  to  interpret  physically  . Another  approach  is  to  compute 
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a one-dimensional  Fourier-  transform  in  terms  of  the  lag,  t,  so  that 
the  resulting  spectrum  is  a function  of  frequency,  f,  and  time,  t. 
This  approach,  however,  has  the  undesirable  effect  that  it  may  be 
negative  at  particular  frequencies  C3 , p . 3G1 3 . 

2.9  Spectral  Estimators 

As  the  spectrum  is  the  Fourier  transform  of  the  autocovariance 
function,  it  is  natural  to  consider  using  the  Fourier  transform  of  a 
sample  autocovariance  function  (such  as  given  in  (2.33)1  as  a 
spectral  estimator;  in  fact,  this  is  commonly  done.  Unfortunately, 
although  sample  autocovariance  functions  are  generally  consistent, 
their  Fourier  transforms  are  not;  the  variance  does  not  tend  to  zero 
for  large  sample  lengths.  Consequently,  smoothing  techniques  have 
been  developed  to  reduce  the  variance  of  spectral  estimators  defined 
in  this  fashion.  Ue  first  describe  the  properties  of  unsmoothed 
spectral  estimators;  smoothed  estimators  are  discussed  in  section 
2.10. 

Assume  x(t)  to  be  a zero-mean,  stationary  random  process  with 
spectrum  G*(f),  and  xT(t)  to  be  a finite  sample  function  of  length 
T.  Then  define  a sample  spectrum  or  periodogram  as 

I„(f)  « J"  (c** 1 » JJc„,,(t) -exp  (-2n  j f t)  dT  (2.40) 
where  ch*{t)  is  the  sample  autocovariance  function  (2.33a)  and 
yiCttyir))  represents  its  finite  Fourier  transform  [4, p. 2151.  As 
discussed  in  Appendix  B,  an  equivalent  definition  of  IK(f)  is 

I*(f)  - (i/t)  |X,(f)  |J  (2.41a) 

and 


xt(f)  » yjx,(tn 


Jo  xT  ( t ) -exp (-2n j f t ) dt 


(2.41b) 


"" 
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(e.g.,  XT(f ) is  the  finite  Fourier  transform  of  xT(t)).  (2.41)  is  a 

particularly  useful  formulation  of  the  periodogram  since  the  Fast 
Fourier  Transform  (see  Appendix  C)  enables  rapid  computation  of 
XT(f)  as  opposed  to  the  relatively  slow  computation  time  for  cxx(t) • 

Ue  are  interested  in  the  expectation  and  variance  of  the 
periodogram.  The  expected  value  may  be  found  by  simply  noting 
E (I* ( f ) ) = j'.;Efcxx(r))-exp(-2njfT)dT 

- J’.;Cxx(t)-(1  - lr|  / T)  *exp  ( — 2n  j f r)  dT  (2.42) 

where  E(cxx(r))  is  given  oy  (2.34a).  Clearly,  because  of  the  finite 
limits  of  integration  and  the  factor  (1  - lr|/T),  the  periodogram  is 
a biased  estimator  of  G>(f).  Note,  however,  that  as  T ■*  ®, 

Via?  E Civ  ( f ) 1 * j;:Cxx(T)-expl-2njfT)dT  = Gx(f)  (2.43) 

so  that  the  periodogram  is  asymptotically  unbiased  and 

E(Ix(f))  * Gx(-:)  . (2.44) 

This  leads  to  an  alternate  definition  of  the  spectral  density 

function  as 

Gx( f ) = Via? E nx(f ) ) = !ia'E((i/T)ixT(f)i*i  . (2.45) 

Ule  can  interpret  (2.42)  as  being  -the -Fourer  transform  of  a 

windowed  autocovariance  function  where 
w,(t)  = (1  - Irl/T).  Irl  < T 

- 0,  otherwise  (2.46). 

is  called  the  Bartlett  window  [5, p.443).  Using  the  fact  that  the 
Fourier  transform  maps  multiplication  into  convolution,  (2.42)  can 
be  written  as  a convolution  of  the  spectrum,  Gx(f),  with  the 

frequency  representat ion  of  (2.4S),  i.e., 

E (Ix(  f ) I = Gx  ( • ) ®U,  (f)  (2.47a) 

and 
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W. ( f ) - f{w,{t)l  - (TMsin(ixTf) /nTf)2  . (2.47b) 

Since  it  is  well  Known  tnat  for  large  T the  representation  of  U,(f) 
in  (2.47b)  tends  to  an  impulse,  we  again  see  that  I x ( f ) is 
asymptotically  unbiased. 

Blackman  and  Tukey  121)  introduced  the  terminology  of  calling  a 
window  in  the  time  domain  (as  in  (2.46))  a lag  window  and  its 
frequency  representation  (2.47b)  a spectral  window.  Several 
different  spectral  windows  have  been  developed  with  various 

properties  (e.g.,  see  (4, p.244)  and  (271).  The  selection  of  an 
appropriate  window  for  a particular  application  is  the  focal  point 
of  much  of  the  research  in  spectral  analysis. 

An  important  consequence  of  (2.42)  and  (2.47)  is  that  the 
smoother  the  spectrum,  G„(f),  the  less  biased  the  estimator  tends  to 
be.  For  a purely  random  process  with  Cxx(r)  = o-2x-S(r),  Ix(f)  is  an 
unbiased  estimator  for  all  T since  [4,p,238] 

E II* ( f ) ) = GxO  - <r2x  . (2.48) 

Conversely,  the  spectrum  tends  to  be  distorted  in  the  vicinity  of 
sharp  peaks  due  to  convolution  with  the . side . lobes  in  the  spectral 
window.  For  this  reasor , it  is  usually  desirable  to  use  a sample 
autocovariance  function  with  its  associated  spectral  window  having 
side  lobes  as  small  as  pcssible.  Alternatively,  bias  may  be  reduced 
if  the  signal  is  pre-whitened  by  passing  it  through  a linear  system 
with  a frequency  response  equal  or  approximately  equal  to  the 

inverse  of  the  anticipated  spectrum  Gx(t). 

General  expressions  for  the  variance  of  the  periodogram  are 
quite  involved  but  have  been  derived  by  Jenkins  and 

Uatts  (4, pp. 412-18)  and  others.  For  a Gaussian,  zero-mean,  white 
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process  x(t)  ~ N(0,<t-2x),  however,  it  can  be  shown  [4, p.2333  that  at 
the  harmonic  frequencies;,  f = k/T,  |kl  = 0,  1,  2,  •••  (and  all 

frequencies  for  large  T) 

var  (Ix( f ) ) = <r\  . (2.43) 

To  a good  approximation  for  non-white  and  non-Gaussian  processes, 
(2.47)  may  be  extended  to  become  [4, p. 2501 

vardy(f))  s G<2(f)*  Cl  + (sin  (2nTf ) / 2itTf ) 2] 

s GxMf)  . (2.50) 

The  exact  formulation  of  (2.50)  is  not  as  important  as  the  fact  that 
it  shows  that  the  variance  does  not  tend  to  zero  for  large  T but 
rather  to  a constant  approximately  equal  to  the  square  of  tn.e 
spectrum  itself.  Ix(f)  is  thus  an  inconsistent  estimator  of  Gx(f). 

In  Appendix  B,  it  is  shown  from  (2.41a)  that  a zero-mean,  white 
Gaussian  process  x(t)  ~ N(0,o-Jx),  the  quantity  2I„(f)/cr2x  is 

distributed  exactly  as  chi-square  with  two  degrees  of  freedom.  For 
other  processes,  if  T is  large,  then  2Ix(f)/Gx(f)  is  approximate  1 y 
X'j.  By  using  the  appropriate  functions  from  Table  2.1,  asymptotic 
results  similar  to  (2.44)  and  (2.50)  are  easily  derived  directly 
from  the  properties  of  the  chi-square  distribution.  In  chapters  3 
and  4 we  make  wide  use  of  the  distribution  of  the  periodogram. 

2.10  Smooth  Spectral  Estimators 

In  194G  Daniel  1 (28)  suggested  that  consistent  estimates  of  the 
spectrum  could  be  obtained  by  averaging  the  periodogram  at  adjacent 
f requenci es . Other  research  by  Bartlett,  Blackman  and  Tukey,  and 

others  [8, p. 2581  extended  and  modified  this  idea  and  helped 
introduce  the  notion  of  a smooth  spectral  estimator. 
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The  Daniel  1 estimator  is  actually  part  of  a broader  class  of 
smooth  spectral  estimators  formed  by  convolving  the  periodogram  with 
a spectral  window  so  that 

S„(f)  ■=  I,(f)©U(f)  = f lcxx(r)-w(t))  . (2.51) 

where  Sx(f)  is  the  smoothed  estimator  and  w(t),  U(f)  are  a 
lag-spectral  window  pair.  As  indicated  in  (2.51),  this  convolution 
is  the  equivalent  of  multiplying  the  autocovariance  estimator  by  an 
appropriate  lag  window.  Uhile  in  theory  any  window  may  be  used,  in 
practice  the  selection  is  usually  limited  to  windows  for  which 
U(f)  > 0 for  all  f.  If  W(f)  is  negative  for  any  f,  Sx(f)  may  also 


be  negative  for  particular  frequencies. 

In  general,  lag  wincows  have  the  properties 

w (0)  = 1.  (2.52a) 

w (-t)  = w ( t) , (2.52b) 

and 

w(t)  =0,  I tl  > H,  H < T . (2.52c) 

Note  especially  that  while  the  sample  length  is  T,  the  window  is 

autocovariance  estimator  need  only  be  computed  for  lags  up  to  M. 

The  expected  value  cf  Sx(f)  is 

E ISX  ( f ) ( = f IEIcxx(r))-w(t)) 


* f ICxx(r)-w,(t)-w(t)i 

= Gx(f)®lJ.(f)®U(f)  2=  Gx(f)  ®U(f ) (2.53) 

where  the  approximation  in  the  last  step  is  for  T much  larger  than  II 
(since  under  that  condition,  U,(f)  will  be  narrow  compared  to  U(f) 
and  much  more  like  an  impulse).  If  the  spectrum  is  sufficiently 
smooth,  then  ElSx(f)l  = C-X(f).  If  G»(f)  is  not  smooth,  then  clearly 
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the  narrower  the  spectral  window,  the  more  fidelity  in  the 
convolution  of  (2.53)  and  the  smaller  the  bias.  In  fact, 
approximate  expressions  of  the  bias  of  Sy; C f ) have  been  computed  for 
various  spectral  windows  [4, p. 2471.  In  general  these  expressions 
are  proportional  to  (1/M)n  where  n > 0 is  some  integer  power. 
Clearly,  then,  as  M increases  (with  a coresponding  decrease  in  the 
bandwidth  of  the  spectral  window),  the  bias  will  be  smaller . From 
this  point  of  view,  then,  it  is  desirable  to  choose  the  width  of  the 
lag  window  to  be  as  large  as  possible. 

Approximate  expressions  of  the  variance  of  Sx(f)  have  also  been 
derived  14, p. 251] , 15, p. 552) 

var  (S„(f)  I ^ (i/T)-G*Mf)©UJ(f) 
s (i/T)-G/(f).j;:UMf)df 

* (i/T)-G/(f)-j;:wMt)dt  * Gx2(f ) • (K  / T ) (2.54) 


where  the  approximation  in  the  last  step  is  again  for  Gx(f)  smooth 
and  N larger  than  M.  This  the  variance  is  now  not  only  proportional 
to  the  square  of  the  spectrum,  but  to  the  area  under  the  squared  lag 


spectral  windows  have  been  computed  (e.g.,  see  [4, p. 2521)  and  are 
generally  proprotional  to  fl.  Thus  we  see  that  the  variance  of  Sx  ( f ) 
is  reduced  as  T increases  and  tf  decreases  (making  Sx(f)  consistent). 
This  last  condition,  however,  is  the  opposite  of  that  required  to 


decrease 

the 

bias. 

As  usual. 

then 

, a 

compromise  must  be 

achieved 

between 

sma  1 1 

variance  (small 

m 

and 

smal  1 

bi3s  (large 

M)  . As 

mentioned  in 

section 

2.3  this 

is  often 

done 

by  minimizing 

the  mean 

scju are  error,  mse(Sx(f)>  = var!S»(f))  + B2i5x(f)l. 

The  above  results  rmy  be  expressed  differently  bu  approximat ing 
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(1)  nS* ( f ) / G* ( f ) with  a chi-square  distribution  having  n = 2T/K 
degrees  of  freedom  (given  by  the  equivalent  degrees  of  freedom 
(2.18a)),  and  (2)  the  effective  bandwidth,  (3,  of  the  spectral  window 
where  (3  * (1/K).  Recall  that  K ~ H and  note  that  n ^ 2T(3  and 

var(Sx(f)l  2!  2G„J(f)/n  ^ G*2(f)/T|3.  As  the  bandwidth  increases, 

then,  the  degrees  of  freedom  becomes  larger  and  the  corresponding 
variance  smaller.  However,  the  increased  bandwidth  means  poorer 

resolution  and  larger  bias.  Conversely,  if  M is  large,  then  the 
bandwidth  is  small  giving  better  resolution  and  smaller  bias,  but 
the  degrees  of  freedom  is  smaller  with  a corresponding  larger 
variance.  The  appropriate  choice  of  a spectral  window  to  compromise 
this  situation  is  a fundamental  topic  of  spectral  analysis 
research  [4, pp. 252-571 . 

In  1348,  Bartlett  (231  proposed  a slightly  different  method  of 
computing  smooth  spectral  estimates.  From  (2.43)  we  see  that  an 
alternate  definition  of  the  spectrum  involves  the  expected  value  of 
the  periodogram.  It  would  seem  logical,  then,  to  improve  the 
periodogram  estimator  by  computing  its  sample  mean.  If  several 
sample  functions  were  available,  this  could  be  done  by  computing 
several  periodograms  and  averaging  frequency  by  frequency;  this  is 
usually  not  the  case,  however.  If  the  process  is  ergodic,  then  by 
applying  the  concepts  of  section  2.G  we  can  compute  per iodograms 
from  adjacent  segments  of  the  one  sample  function  xT(t).  These 
peridograms  are  then  averaged  to  produce  Bartlett’s  smooth 
estimator  (4 , pp. 233-431 . Thus  we  define 

P*(f)  - l/dl."j.(f)  (2.55a) 


and 
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Ii.( f ) = (i/T)  • I J”'xi(t)-exp(-2njft)cltlI  (2.55b) 

where  x(t)  has  been  divided  into  N segments  of  length  II  such  that 
M*N  = T and  I;(f)  is  the  periodogram  computed  from  the  ith  segment, 
xL(t)  .+ 

The  mean  and  variance  of  Px(f)  come  directly  from  (2.47), 
(2.50),  and  (2.55a).  If  x(t)  is  Gaussian  and  white,  then  x,(t)  is 

independent  of  x;(t),  i * j and,  hence,  I;(f)  is  independent  of 
I j ( f ) . If  the  process  is  correlated,  but  Cxx(r)  is  small  for  r > H, 
the  segment  length,  then  It(f)  is  still  nearly  independent  of  Ij(f). 
Thus  it  is  reasonable  to  assume  that  cov  (I, ( f ) , Ij ( f ) ! =0,  i * j , so 
that 

E )PX(  f ) 1 = i/NliNiEII,(f))  = Eil.(f))  = G„(f)®U,(f)  (2.56a) 

and 

variPx(f))  = ( i/N) 2,Nivar  tl,(f)l 

= varlljf))  /N  s Gx2(f) /N  (2.56b) 

where  U,(f)  is  the  BaMlett  window  (2.47b)  on  [-H,  Ml.  For  a 
Gaussian,  zero-mean,  while  process  x ( t ) ~ N(0,<72x),  (2.56)  becomes 

E (P«(f ) t = a2,  (2.57a) 

and 

var  )PX ( f ) ) = tr\/N  . (2.57b) 

Like  Ix(f),  Px(f)  is  asymptotically  unbiased.  However,  it  is  a 
consistent  estimator  since  varlP„(t)l  tends  to  zero  for  large  N. 

The  tradeoff  between  bias  and  variance  encountered  with  the 


+Note  that  the  frequency  smoothing  procedure  of  Daniel  1 also 
computes  a form  of  the  sample  mean  of  Ix(f)  and  is  often  called 
frequency  smoothing.  Bartlett’s  procedure  is  similarly  called  time 
or  ensemble  smoothing. 
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previously  discussed  class  of  smoothed  estimators  also  apply  to 
PK ( f ) . As  before,  an  effective  bandwidth  and  equivalent  degrees  of 
freedom  can  be  computed.  From  the  relationship  n a;  2T(3  where  n is 
the  equivalent  degrees  cf  freedom  and  (3  is  the  bandwidth,  we  have 


(3  a*  n / 2T. 

I n 

Appendix 

A, 

however,  it  is 

shown  that  n is 

approximately 

2N 

(exac  1 1 y 

2N 

for  a Gaussian, 

white,  zero-mean 

process)  and  thus  (3  s 2N/2T  = 1/H.  For  finite  data,  T is  fixed. 
Thus,  as  M is  made  smaller,  N becomes  larger  and  the  variance  is 
reduced.  However,  this  again  produces  less  resolution  with  a 
subsequent  increase  in  bias.  Conversely,  fewer  periodograms 
improves  the  resolution  and  bias,  but  increases  the  variance. 

The  Bartlett  procedure  is  particularly  useful  in  practice  since 
it  enables  periodograms  to  be  computed  for  short  segments  of  data. 
Ulith  application  of  the  Fast  Fourier  Transform  (see  Appendix  C)  , 
these  computations  mat  be  made  rapidly  and  the  resulting 
periodograms  averaged.  For  this  reason,  the  remainder  of  this 
effort  centers  on  the  Bartlett  estimator  with  some  slight  but 
significant  modification;  to  produce  our  log  spectral  estimators. 

A modified  form  of  the  Bartlett  procedure  was  proposed  by 
Uelch  [221  in  1967.  In  this  case,  a lag  window  is  applied  to  each 
of  the  data  segments,  x^t),  before  computing  the  Fourier  transform 
giving 

P„(f)  - l/tiLV.ff),  (2.58a) 

J;(f)  = (1/ m.J)].|J"i;xl(t)-w(t)-exp(-2njft)dt|J,  (2.58b) 

and 

U = (1  /tl)  (2.58c) 

where  the  factor  U insures  that  J„  ( f ) is  asymptotically  unbiased. 
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Welch  gives  the  expected  value  and  variance  as 

E (Px(f)  1 = Gx(;)®W(f)  , (2.59a) 

U(f)  = Cl  / (n-U)  3 *|  J'oW  ( t)  *exp  (~2n  j f t)  dt  I2,  (2.59b) 

and 

van  lPx(f) ) * G,2(f)  / N (2.59c) 

where  U ( f ) is  now  the  magnitude  squared  of  the  Fourier  transform  of 
the  lag  window. 

This  procedure  enables  selection  of  an  appropriate  spectral 
window  with  desired  properties,  yet  retains  the  computational 
advantages  of  the  Bartlett  estimator.  Note  in  particular  that 
windows  with  negative  values  may  be  used  since  the  expected  value 
(2.59a)  now  involves  convolution  with  the  square  of  the  magnitude  of 
the  spectral  window.  Welch  extends  his  analysis  to  the  case  where 
the  data  segments,  x,(t),  overlap.  Although  the  periodograms  are  no 
longer  independent,  thus  increasing  the  variance,  more  are  available 
to  be  averaged.  Welch  fcund  that  this  is  a net  gain  in  reduction  of 
the  variance.  In  this  research,  we  are  concerned  only  with 

non-overlapping  segments  and  the  corresponding  assumption  that  the 
periodograms  are  independent. 

Welch’s  algorithm  is  widely  used  in  modern  spectral  analysis 

and  is  the  one  used  far  the  experimental  computations  in  this 

research.  Further  details  of  the  digital  implementation  of  (2.58) 
are  given  in  Appendix  C including  a discussion  of  the  spectral 

window  used. 

Figures  2.2  and  2.3  illustrate  the  effect  of  spectral 
smoothing.  Figure  2.2(a)  is  a periodogram  computed  from  simulated 
Gaussian,  zero-mean  white  noise  with  <r\  = 1000  (see  Appendix  C for 


a detailed  description  of  the  production  of  this  data).  It  is 
apparent  from  this  graph  that  ElPy(f)l  = EII*(f))  = 1000  and 

varllx(f)l  ■=  (1000)2  as  predicted  by  (2.57).  Figure  2.2(b)  is  the 
average  of  4 periodograms  for  this  same  data.  Note  how  the  standard 
deviation  has  been  reducsd  by  a factor  of  2 (since  the  variance  is 
reduced  by  1/N,  the  standard  deviation  is  reduced  by  (1  / N) i/z)  . 
Similarly,  figures  2.2(ct  and  2.2(d)  arc  smooth  estimates  with  N = 
IB  and  G4  respectively.1- 

Figure  2.31-1-  parallels  figure  2.2  for  a colored  Gaussian 
process  (the  white  process  depicteci  in  figure  2.3  was  filtered  using 
the  linear  system  shown  in  figure  4.8(c)).  Note  how  for  N = 1 

(figure  2.3(a))  the  stancard  deviation  is  as  large  as  the  signal  and 
masks  all  but  the  general  shape  of  the  spectrum.  After  smoothing  B4 
periodograms  (figure  2.3(d)),  the  nature  of  the  spectrum  is  much 
more  evident. 


+The  standard  deviation  of  this  and  other  figures  may  be  visually 
judged  as  being  half  the  "fuzziness"  surrounding  the  apparent 
expected  value.  The  variance  is  then  the  square  of  the  standard 
deviation. 

++In  figure  2.3  and  the  remaining  figures,  unless  otherwise  noted, 
the  abcissa  represents  frequency  in  Hertz  on  a logarithmic  scale. 
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FIGURE  2.3 

Progressively  smoothed  spectral  estimates  for  a Gaussian,  colored 
process,  11(0,1000)  with:  (a)  N = 1 (periodogram) , (b)  N = 4,  (c)  N = 
16,  and  (d)  N = 64. 
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CHAPTER  3 


LOG  SPECTRAL  ESTIMATES  FOR  STATIONARY  PROCESSES 

3.1  Logarithmic  Representation  of  Data 

Physical  data  are  commonly  displayed  on  a logarithmic  scale. 
This  is  because  (1)  the  logarithm  compresses  the  dynamic  range 
permitting  small  changes  in  the  value  of  the  data  to  be  observed 
simultaneously  with  large  ones  as  well  as  enabling  digital  storage 
with  fewer  bits;  (2)  such  a transformation  may  produce  data  which  is 
more  nearly  Gaussian  in  character  [30, p. 741;  and  (3)  as  will  be 
shown,  it  can  result  in  c.ata  with  a variance  or  confidence  intervals 
independent  of  the  value  of  the  signal.  Examples  include  the 

logarithmic  Richter  scale  used  in  seismology  and,  of  course,  power 
spectra  represented  in  decibels. + 

3.2  Log  Spectral  Density  Function 

In  section  2.8  the  spectral  density  function  was  defined  for  a 
stationary  random  process.  I t is  a simple  extension  of  this  concept 
to  define  the  log  spectral  density  function  (log  spectrum)  as 


'’"Log  spectra  are  commonly  presented  on  the  decibel  scale  (dB  » 
10> log,0G  where  G represents  power).  Throughout  these  discussions, 
however,  derivations  are  in  terms  of  the  natural  logarithm  to 
produce  equations  of  manageable  mathematical  complexity.  Data  and 
figures,  though,  are  in  cecibels  unless  othewise  noted.  In  all 
cases,  results  of  the  equations  can  be  converted  to  decibels  by 
using  the  appropriate  multiplicative  constant  derived  from  the 
identity  logl0x  - logex/log,10  - (0.43429-”) logex. 
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e*m  - logGx(f ) » logf{Cxx(r))  (3.1) 

where  Cxx(t)  i3  the  autocovariance  function  associated  with  the 
stationary  process  x(t).  Since  the  spectrum  is  non-negative,  the 
log  spectrum  is  defined  for  all  frequencies  except  where  Gx(f)  » 0. 

3.3  Log  Spectral  Estimation 

In  sections  2.9  and  2.10  we  derived  formulas  for  the 
expectation  and  variance  of  various  spectral  estimators.  He  now 
extend  these  concepts  oy  defining  log  spectral  estimators  and 
describing  their  statistical  properties. 

An  obvious  log  spectral  estimator  is  simply  the  logarithm  of 
the  periodogram,  1 1 ( f ) . Similarly,  we  could  compute  the  logarithm 
of  one  of  the  smoothed  spectral  estimators.  Alternatively,  we  could 
effect  the  logarithmic  transformation  before  the  smoothing  process. 
As  we  shall  see,  reversing  the  ordering  of  the  logarithmic 
transformat  ion  and  smoothing  process  produces  two  significantly 
different  estimators  (particularly  for  nonstationary  processes). 

This  research  is  limited  to  a discussion  of  two  particular  log 
spectral  estimators.  Both  are  smoothed  by  the  Bartlett  averaging 
procedure  of  non-overlapping  periodograms  or  their  logarithms.  Log 
spectral  estimators  smoothed  in  frequency  by  convolution  with  a 
spectral  window  are  not  explicitly  considered. 

3. A Log  Average  Spectrum 

The  log  average  spectrum  (LAS)  is  the  logarithm  of  the  smoothed 
spectral  estimator,  P„(f)  (2.55a).  He  thus  have 

px(f)  * logP.(f)  *>  log  C (UN)  2.*ili  ( f ) ) (3.2) 


where  P„(f)  is  the  log  average  spectrum,  and  I,(f)  is  given  by 
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(2.55b) . 

3.5  Average  Log  Spectrum 

The  average  log  spectrum  (ALS)  is  derived  by  computing  the 
logarithmic  average  of  the  periodograms.  Ue  thus  have 

L*( f ) - l un)  = li/NJZAlofllitf)  (3.3) 

where  Lx(f)  is  the  average  log  spectrum. 

As  is  discussed  more  fully  in  chapter  5,  the  motivation  for 
this  estimator  conies  from  consideration  of  a generalized  linear 
system  that  obeys  superposition  across  convolution  [161 . Such  a 

system  uses  a Fourier  transform  followed  by  a logarithm  to  map 
convolution  into  addition.  At  this  point,  linear  operations  may  be 
performed.  If  this  operation  is  an  averaging  process,  and  only 
magnitudes  are  considered,  this  procedure  is  equivalent  to  computing 
the  average  log  spectrum. 

3.6  Statistical  Properties  of  the  Log  Average  Spectrum 

The  difference  between  linear  and  logarithmic  averaging  has 
been  discussed  in  the  literature  for  a variety  of  distributions  and 
applications.  Cox  [91  shows  that  this  difference  has  a lower  bound 
which  is  a function  only  of  the  dynanmic  range  of  the  data  (this  is 
discussed  further  in  section  3.9).  His  result  is  independent  of  any 
soecific  distribution.  Mitchell  (11)  has  derived  expressions  for 
the  expectation  of  this  difference  for  four  distributions,  including 
the  uniform  and  lognormal;  Hershey  (10)  for  the  expectation  and 
variance  assuming  Gaussian  data;  and  Musal  (12)  and  Sugai  and 
Christopher  [14]  for  the  Rayleigh  and  Maxwell  distributions.  In  a 
recent  paper,  Ricker  and  Uilliams  [13]  discuss  the  advantages  of  log 


40 


power  estimates  in  terms  of  the  logarithmic  average  of  chi-square 
data.  Results  similar  to  those  of  Ricker  and  Uilliams  form  the 
basis  for  the  remainder  of  this  work. 

The  statistical  properties  of  (3.2)  and  (3.3)  are  most  easily 
derived  if  x(t)  is  a stationary,  Gaussian,  white  process  with 
spectral  density  function  Gx(f)  = cr2x.  Departures  from  these 
conditions  are  discussed  in  section  3.9. 

With  these  assumptions,  as  discussed  in  Appendix  B, 
2N.P* ( f ) / crJx  is  distributed  as  x2,„  (for  N averaged  periodograms)  and 
thus  px(f)  « logP><(f)  = log  (^Vx^/  2N ) . From  Table  2.1,  with  r ■= 
tr%/2N  and  n = 2N,  we  have 

E(px(f))  - logcr2*  + i|HN)  - log(N)  a log(7-2x  (3.4a) 

and 

var  (px ( f ) ) = 4j‘  (N)  a 1/N  (3.4b) 

where  4/(N)  is  the  digamma  function  (see  Appendix  A).  For  large  N 
(e.g.,  N > 20)  4/  (N)  a loy(N),  4/'  (N)  a 1/N  and  the  approximations  of 
(3.4)  hold. 


3.7  Statistical  Properties  of  the  Average  Log  Spectrum 

Under  the  assumptions  of  section  3.G,  2It(f)  / o-2*  is  distributed 
as  X%  and  so  ix(f)  ■ logl^f)  ■*  log(crVx??/2) . Again,  using  Table 
2.1,  with  r = <r\!2  and  n - 2,  and  noting  that  IJf)  is  independent 
of  I; ( f ) , i * j 

E !Lx(f) ) = E ( (i/N)  2,NtlogI,(f)  I = (i/n)  I£E  (log  Ij.  ( f ) ) 

« (i/n)  ZAlog^  - y - logtr2,  - y (3.5a) 

and 

var(Lx(f)l  = var  t (i/n)  lAlogl^fH 


(l/N) *2i\var  tloglitf ) » = nV  (BN) 
n1/  (8N)  = 1.6449—/N 
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(3.5b) 

where  v is  Euler’s  constant  ( V = 0.57721*..). 

From  (3.4)  and  (3.5)  we  can  see  that  for  stationary,  white 
Gaussian  processes,  the  two  estimators  are  essentially  equivalent 
since  their  expected  values  differ  (asymptotically)  by  a universal 
constant  independent  of  the  process.  The  log  average  spectrum  is 
asymptotically  unbiased;  the  average  log  spectrum  biased.  However, 
this  latter  bias,  being  constant,  is  easily  removed. 

For  both  estimators,  var(px(f)i,  var(Lx(f))  -*  0 as  N -»  «. 
However,  for  a given  value  of  N,  clearly  var(Lx(f))  > varlpx(f)). 
Specif ical ly , 

var (Lx ( f ) ) / var  ipx ( f ) 1 

= (itJ/6)/^'(N)  s 1.6449... /N  . (3.6) 

Thus  the  variance  of  the  log  average  spectrum  is  about  39%  less  than 
that  of  the  average  log  spectrum.  Equivalently,  the  log  average 
spectrum  has  22%  less  standard  deviation  (is  22%  smoother) . 

3.8  The  Activation  Spectrum 

A useful  and  interesting  finding  of  this  research  is  the 
difference  between  the  ALS  and  LAS  estimators 

Ax(f)  = Ly(f)  -•  px(f)  (3.7) 

where  Ax(f)  is  termed  the  activation  spectrum1".  From  (3.4)  and 
(3.5)  we  see  that  the  expected  value  of  Ax(f)  is  given  by 


+The  term  activation  spectrum  was  suggested  to  the  author  by  T. 
G.  Stockham,  Jr.,  Department  of  Computer  Science,  University  of 
Utah.  Motivation  for  this  term  is  discussed  in  section  4.4. 


E (A*  ( f ) ) = log  (N)  - *1/  (N) 
For  a stationary  process,  then, 


y a*  -y  , 
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(3.8) 

the  activation  spectrum  is 
aymptotical ly  equal  to  Euler’s  constant. 

As  discussed  in  chapter  4,  (3.8)  does  not  hold  for 
nonstationary  processes.  In  fact,  the  activation  spectrum  may  have 
a value  significantly  different  from  Euler’s  constant,  providing  a 
test  for  nonstationari ty  in  the  presence  of  stationary  noise. 

Even  for  a stationary  process,  (3.8)  is  particularly  useful. 
For  example,  if  it  is  desired  to  identify  a coherent  signal  in 
stationary  backgrond  noise,  there  will  be  a higher  coherent  signal 
to  noise  ratio  in  the  ALS  estimator.  Since  the  coherent  peak  will 
have  nearly  the  same  val  ue  in  both  estimators  (not  being  random)  , 
the  background  noise  will  be  down  an  additional  2.5  dB  (the  decibel 
equivalent  of  Euler’s  constant)  in  the  ALS  estimator.  However,  this 
increase  in  S/N  ratio  is  achieved  at  the  expense  of  stability  (i.e., 
greater  variance) . 

For  a given  set  of  periodograms,  the  activation  spectrum  is  the 
logarithm  of  the  ratio  of  their  geometric  to  arithmetic  means.  It 
is  well  known  that  for  any  data,  this  ratio  has  an  upper  bound  of  1 
and,  thus,  its  logarithm  has  an  upper  bound  of  0.  In  19),  Cox  shows 
that  this  ratio  has  e lower  bound  which  is  a monotonical  1 y 
decreasing  function  (only)  of  the  dynamic  range  or  activation  of  the 
data.  Specifically,  if  X = max  IIk(f) ) /min  (I  t ( f ) ) is  the  dynamic 
range  of  a set  of  periodograms  (termed  the  spectral  dynamic  range), 
then 

0 > A,(f)  > logB(K)  (3.9) 


where  logB(K)  is  given  by 
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log  EB  (K) ] - 1 + log  [log  (K) / (K  - 1)1 

- log (K) / (K  - 1)  . (3.10) 

This  relationship  is  usee,  again  in  chapter  4. 

3.9  Results  For  Non-Gau= sian,  Correlated  Processes. 

The  results  of  the  previous  sections  are  exact  for  stationary, 
Gaussian,  zero-mean,  white  processes.  In  practice,  these  conditions 
are  seldom  if  ever  met  and  the  effects  of  departures  are  important. 
In  chapter  4 we  discuss  nonstationarity  while  in  this  section, 
non-Gaussian,  correlated  processes  are  considered. 

The  restriction  to  2 zero-mean  process  is  not  severe.  It  may 

be  easily  met  by  subtracting  the  sample  mean  of  the  process.  If  not 

clone,  the  resulting  spectrum  will  have  an  impulse  at  f = 0. 

Resulting  spectral  and  log  spectral  estimators  may  then  tend  to  be 
heavily  biased  in  the  vicinity  of  the  origin  due  to  convolution  with 
side  lobes  of  the  spectre, 1 window. 

The  restriction  to  a Gaussian  process  is  similarly  not  severe. 
If  the  process  is  not  highly  correlated,  its  Fourier  transform  will 
consist  of  sums  of  essentially  independent  random  variables.  By  the. 
Central  Limit  Theorem,  then,  X,(f)  ° flx^tt))  will  tend  to  normality 
regardless  of  the  distribution  of  x,(t);  this  is  particularly  true 
for  large  segment  lengths,  M. 

Correlation  in  the  process,  however,  may  be  more  significant. 
If  the  process  is  highly  correlated,  the  resulting  per iodogranis  will 
be  biased  in  accordance  with  (2.47-3) . E!I;(f))  will  depart  from 

Gy(f)  as  a function  of  the  spectral  window,  U8(f),  and  the  shape  of 
the  spectrum.  This  may  be  thought  of  as  follows.  In  the  vicinity 
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of  peaks  in  the  spectrum,  convolution  with  the  spectral  window  has 
the  effect  of  increasing  the  degrees  of  freedom  of  the  chi-square 
random  variable  associated  with  each  frequency.  Thus,  they  may  no 
longer  be  distributed  as  x??,  but  as  x2„  with  n given  by  the 

equivalent  degrees  of  freedom  (2.18a).  Conversely,  X„(f)  and  X,(f) 
(the  real  and  imaginary  parts  of  X^f),  respectively)  may  no  longer 
be  independent,  thus  reducing  the  degrees  of  freedom.  Clearly, 
these  considerations  will  affect  (3.4)  and  (3.5).  Also,  if  the 
process  is  highly  correlated,  the  Ii(f)s  may  not  be  independent  and 
their  covariance  would  be  a significant  term  in  (3.5b).  A good 

example  of  this  latter  point  is  discussed  by  Welch  (221  in  which  he 

considers  overlapping  periodograms  that  are  definitely  not 

independent. 

I f the  spectrum  is  reasonably  smooth  and  N and  T are  large, 
then  the  results  for  Gaussian,  white  processes  may  be  extended  with 
reasonable  accuracy  to  more  general  processes.  In  this  case,  (3.4) 
becomes 

E (Px  ( f ) ) * lcgG*(f ) + i/  (N)  - log  (N)  * logG*(f)  (3.11a) 

and 

var  Ip* ( f ) I * *'(N)  ^ 1/N  (3.11b) 

and  (3.5)  becomes 

EILx(f))  a looG,(f)  - y (3.12a) 

and 

var  (Lx ( f ) I at  n!/  (GN)  * 1.S449-../N  (3.12b) 

where  G*(f)  is  the  spectrum  of  x(t).  Similarly,  for  the  activation 
spectrum,  (3.8)  becomes 


E IAx(f) ) as  leg  (N)  - 4/  (N) 


Y 


-Y  . 


(3. 13) 
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The  basic  relationships  between  the  LAS  and  ALS  estimators  is 

basically  preserved  even  for  more  general  processes. 

3.10  Experimental  Results 

Figure  3.1  illustrates  both  the  log  average  spectrum  (a)  and 
average  log  spectrum  (b)  for  a computer  simulated,  stationary, 
white,  Gaussian  process  x(t)  ~ N(0,1000)  (see  Appendix  C) . From 
(3.4)  and  (3.5)  we  would  expect  p*(f)  = 30  dB,  and  L*(f)  = 27.5  dB. 
These  values,  as  well  as  the  greater  stability  of  the  LAS  estimate 
are  clearly  observed.  Figure  3.1(c)  illustrates  the  activation 
spectrum  and  figure  3.1(d)  the  spectral  dynamic  range  for  this  same 
data.  The  predicted  value  of  the  activation  spectrum,  Ax(f)  ■=  -2.5 
dB,  i9  apparent.  These  estimates  were  computed  by  averaging  200 

periodograms. 

Figure  3.2  is  similar  to  figure  3.1  but  is  for  a colored, 
stationary  process.  In  this  case,  the  stationary,  Gaussian  data  of 
Figure  3.1  has  been  passed  through  a system  with  a non-flat 
frequency  response  (see  figure  4.8(c))  and  the  log  spectrum 

estimated.  As  before,  -igure  3.2(a)  is  the  LAS  estimator,  figure 
3.2(b)  the  ALS,  figure  3.2(c)  the  activation  spectrum  and  figure 
3.2(d)  the  spectral  dynan ic  range.  Note  that  although  this  is  now  a 
correlated  process,  the  activation  spectrum  is  still  centered  about 
-2.5  dB  as  predicted  by  (3.13). 

Tables  3.1  - 3.4  present  the  sample  mean  (2.22a)  and  variance 
(2.22b)  for  the  data  of  figures  3.1(a)  and  3.1(b)  (for  various 
numbers  of  averaged  per lodograms) . Also  given  are  the  predicted 

values  from  (3.4)  and  (3.5)  as  well  as  well  as  confidence  intervals 
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05%)  computed  using  (2.23)  and  (2.24).  Note  again  that  these 
figures  are  decibel  equivalents,  e.g.,  V - 2.5  dB.  By  comparing 
values  it  is  clear  that  there  is  good  agreement  between  experiment 

F'  / 

and  theory. 


■*-%t  **>*•%■• 


FIGURE  3.1 

Log  spectral  estimates  for  a Gaussian,  white  process,  N(0,1000) 
(linear  x-axis):  (a)  log  average  spectrum  (N  = 200),  (b)  average  log 
spectrum,  (c)  activation  spectrum,  and  (d)  spectral  dynamic  range. 


FIGURE  3.2 


i*\  i*~°^  spectra^  estimates  for  a Gaussian,  colored  process,  N ( 0 , 1 000 ) • 
a log  average  spectrum  (N  = 470),  (b)  average  log  spectrum  (N  = 470) 
(c)  activation  spectrum,  and  (d)  spectral  dynamic  range. 
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TABLE  3.1 

SAMPLE  MEAN  CF  THE  LOG  AVERAGE  SPECTRUM  FOR  A 
STATIONARY.  GAUSSIAN.  WHITE  PROCESS  ~ N (0,1000) 


N 

Loner 
Conf id. 

Sample 

Mean 

Upper 
Conf id. 

Predict 

Value 

1 

26.5937 

27.5184 

28.4431 

27.4932 

2 

28.7438 

28.8484 

28.9530 

28.8258 

4 

29.3746 

29.4460 

29.5174 

29.4346 

S 

29.5582 

29.6161 

29.6740 

29.6281 

8 

29.6787 

29.7275 

29.7763 

29.7229 

10 

29.7410 

29.7848 

29.8286 

29.7792 

20 

29.8810 

29.9112 

29.9414 

29.8905 

30 

29.9130 

29.9373 

29.9616 

29.9272 

40 

29.9277 

29.9486 

29.9695 

29.9455 

50 

29.9487 

29.9674 

29.9861 

29.9564 

100 

29.9655 

29.9790 

29.9925 

29.9782 

200 

29.9848 

29.9943 

30.0038 

29.9891 

400 

29.9936 

30.0003 

30.0070 

29.9946 

TABLE  3.2 


SAMPLE  MEAN 

CF  THE  AVERAGE 

LOG  SPECTRUM 

FOR  A 

STATIONARY, 

GAUSSIAN,  WHITE 

PROCESS  ~ N(0 

,1000) 

N 

Lower 

Sample 

Upper 

Predict' 

Conf id. 

Mean 

Conf id. 

Value 

1 

26.5937 

27.5184 

28.4431 

27.4932 

2 

27.4351 

27.5540 

27.6729 

27.4932 

4 

27.4501 

27.5379 

27.6197 

27.4932 

6 

27.4242 

27.4949 

27.5656 

27.4932 

8 

27.4681 

27.5237 

27.5893 

27.4932 

10 

27.4540 

27.5090 

27.5640 

27.4932 

20 

27.4780 

27.5166 

27.5552 

27.4932 

30 

27.4741 

27.5049 

27.5357 

27.4932 

40 

27.4671 

27.4939 

27.5207 

27.4932 

50 

27.4814 

27.5051 

27.5288 

27.4932 

100 

27.4815 

27.4986 

27.5157 

27.4932 

200 

27.4902 

27.5024 

27.5146 

27.4932 

400 

27.4926 

27.5013 

27.5100 

27.4932 

50 


TABLE  3.3 

SAMPLE  VARIANCE  OF  THE  LOG  AVERAGE  SPECTRUM  FOR  A 
STATIONARY,  GAUSSIAN,  WHITE  PROCESS  ~ N (0,1000) 


N 

Lower 
Conf id. 

Sampl e 
Mean 

Upper 
Conf id. 

Predicti 

Value 

1 

28.9348 

30.2015 

31.5535 

31.0254 

2 

11.1743 

11.6635 

12.1857 

12.1642 

4 

5.1990 

5.4266 

5.6696 

5.3532 

6 

3.4275 

3.5775 

3.7377 

3.4200 

8 

2.43G8 

2.5434 

2.6573 

2.5111 

10 

1 . 9578 

2.0435 

2.1350 

1.9836 

20 

0.9345 

0.9754 

1.0191 

0.9670 

30 

0.6006 

0.6263 

0.6550 

0.6393 

40 

0.5561 

0.4656 

0.4864 

0.4775 

50 

0.3555 

0.3711 

0.3877 

0.3810 

100 

0.1871 

0.1952 

0.2040 

0. 1896 

200 

0.0913 

0.0953 

0.0996 

0.0945 

400 

0.0457 

0.0477 

0.0498 

0.0472 

TABLE  3.4 


SAMPLE  VARIANCE  OF  THE  AVERAGE  LOG  SPECTRUM  FOR  A 
STATIONARY,  GAUSSIAN,  WHITE  PROCESS  ~ N (0,1000) 


N 

Lower 
Conf id. 

Sample 

Mean 

Upper 
Conf id. 

Predict^ 

Value 

1 

28.9348 

30.2015 

31.5535 

31.0254 

2 

14.4477 

15.0802 

15.7553 

15.5127 

4 

7.3385 

7.6597 

8.0026 

7.7564 

6 

5.1030 

5.3264 

5.5649 

5.1709 

8 

3.7444 

3.9084 

4.0833 

3.8782 

10 

3.0905 

3.2258 

3.3702 

3. 1025 

20 

1.5178 

1.5842 

1.6552 

1.5513 

30 

0.9702 

1.0127 

1 . 0580 

1.0342 

40 

0.7330 

0.7651 

0.7994 

0.7756 

50 

0.5750 

0.6002 

0.6270 

0.6205 

100 

0.2991 

0.3122 

0.3262 

0.3103 

200 

0.1509 

0.1575 

0.1646 

0.1551 

400 

0.0772 

0.0806 

0.0842 

0.0776 

CHAPTER  4 


LOG  SPECTRAL  ESTIMATES  FOR  NONSTATIONARY  PROCESSES 

4.1  The  Problem  of  Nonstationarity 

When  the  tog  spectral  estimates  described  in  chapter  3 are 

computed  for  nonstationary  processes,  results  uith  properties  quite 
different  from  (3.11)  and  13.12)  are  observed.  Since  most  practical 
signals  exhibit  some  nonstationarity.  it  is  important  to  understand 
the  reason,  for  these  differences  and  to  try  to  quantify  them. 

An  inspection  of  figures  4.1  - 4.5  clearly  reveals  some  of 

these  differences. 4 Figure  4.1  presents  estimates  of  the  log 
spectrum  of  a 1907  recording  of  Enrico  Caruso  singing  Vesti  la 
Giubba. " Bofb  the  log  average  spectrum  (a)  and  .be  average  log 
spectrum  (bl  are  eboun  along  uith  the  activation  spectrum  (c)  and 
the  spectral  dynamic  range  (d) . Figures  4.2  and  4.3  are 
estimates  of  mors  recent  recordings  of  the  same  selection  sung  by 
Jussi  Bjoerling  and  Mario  Lanpa.  respectively.  Figures  4.4  and  4.S 
are  log  spectral  sstima.es  for  a female  singer  (Sharon  BrocKbank) 
and  a 3tring  ensemble  digitised  directly  during  live  rec  9 

sessions. 

Tuo  initial  observations  are  apparent.  First,  .be  log  average 
spectrum  has  a considerably  different  shape  from  the  average  log 


tin  these  figures,  the  log  spectral  est 
biased  by  -90.30899872  dE  before  displaying. 


timates  have  been  arbitarily 


FIGURE  4.1 


Log  spectral  estimates  (N  = 470)  for  a 1907  recording  of  Enrico 
Caruso  singing  "Vesti  la  Giubba":  (a)  log  average  spectrum,  (b)  average 
log  spectrum,  (c)  activation  spectrum,  and  (d)  spectral  dynamic  range. 
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FIGURE  4.2 

Log  spectral  estimates  (N  = 544)  for  a modern  recording  of  Jussi 
Bjoerling  singing  "Vesti  la  Giubba":  (a)  log  average  spectrum,  (b) 
averaqe  loq  spectrum,  (c)  activation  spectrum,  and  (d)  spectral  dynamic 
range. 
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FIGURE  4.3 

Log  spectral  estimates  (N  = 390)  for  a modern  recording  of  Mario 
Lanza  singing  "Vesti  la  Giubba":  (a)  log  average  spectrum,  (b)  average 
log  spectrum,  (c)  activation  spectrum,  and  (d)  spectral  dynamic  range. 


e*ea 


1 00or ♦ i 


FIGURE  4.4 


Log  spectral  estimates  (N  = 500)  for  a female  singing  (digitized 
live):  (a)  log  average  spectrum,  (b)  average  log  spectrum,  (c)  activa 
tion  spectrum,  and  (d)  spectral  dynamic  range. 
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FIGURE  4.5 

Log  spectral  estimates  (N  = 50G)  for  a string  ensemble  (digitized 
live):  (a)  log  average  spectrum,  (b)  average  log  spectrum,  (c)  activa- 
tion spectrum,  and  (d)  spectral  dynamic  range. 
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spectrum  for  the  same  process.  Their  difference  is  clearly  greater 
than  2.5  dB.  Secondly,  while  one  would  expect  the  log  average 
spectrum  to  be  more  stable,  in  fact  it  appears  that  the  average  log 
spectrum  is  smoother.  In  the  remainder  of  this  chapter,  we  propose 
a model  for  nonstationary  processes  and  subsequently  explain  these 
observations. 

Before  proceeding,  however,  an  important  point  must  be  clear. 
The  Bartlett  spectral  estimator  developed  in  section  2.10  and  used 
in  chapter  3 estimates  the  single  spectrum  associated  with  a 
stationary  process.  Further,  it  explicitly  assumes  the  process  to 
be  ergodic  so  that  adjacent  data  segments  can  be  used  to  represent 
sample  functions  of  the  ensemble.  Clearly,  both  of  these 

requirements  fail  for  nonstationary  processes.  Not  only  is  the 
process  not  ergodic,  but  because  of  the  nonstationarity,  there  is  no 
single,  mathematically  defined  spectrum  with  which  to  compare  the 
resulting  spectral  estimates.  Hence,  in  these  discussions  we  are 
not  able  to  do  more  than  try  to  understand  what  the  log  average 
spectrum  and  average  log  spectrum  represent  as  computational 
procedures  and  how  to  interpret  them.  It  is  meaningless  to  ask 
which  is  a better  estimator  in  the  sense  that  they  are  not  really 
estimating  a conventional  statistical  parameter  of  the  random 
process. 

Nonetheless,  computation  of  the  average  log  spectrum  and  log 
average  spectrum  can  be  cuite  useful  for  practical  signals.  In  some 
sense,  they  represent  an  "average'  description  of  the  distribution 
of  power  with  frequency  and,  as  such,  estimate  this  "average" 
spectrum.  The  LAS  and  ALS  estimators  differ  in  how  they  compute 


J 


58 


this  average.  In  the  LAS,  small  values  of  the  periodograms 
contribute  little  to  the  average,  while  in  the  ALS,  the  logarithmic 
transformation  results  in  both  large  and  small  values  contributing 
more  equally.  By  properly  interpreting  this  average,  a general 
characterization  of  the  process  is  possible. 

There  is  another  interesting  application  of  these  estimators 
which,  in  fact,  provides  a method  of  directly  comparing  them.  If 
the  nonstationary  process  has  been  passed  through  a linear  system, 
it  is  possible  to  estimate  the  frequency  response  of  this  system 
from  the  log  average  or  average  log  spectra  of  the  input  and  output 
processes.  Equation  (11.38)  shows  the  relationship  between  the 
magnitude  of  the  frequency  response  of  the  linear  system,  H(f),  and 
the  spectra  of  the  input  and  output.  If  we  assume  that  this 
relationship  is  approximately  true  for  a nonstationary  process  with 
the  spectrum  being  replaced  by  the  "average"  spectrum,  then  dividing 
the  "average"  spectral  estimate  of  the  output  by  that  of  the  input 
yields  an  estimate  of  IH(f)l2.  A1  ternatively,  the  ALS  or  LAS 
estimators  may  be  subtracted  to  yield  an  estimate  of  log|H(f)|2. 
In  this  case,  it  is  meaningful  to  ask  which  yields  the  most  stable 
and  the  least  biased  estimate.  This  question  is  explored  more  fully 
in  chapter  5. 

A. 2 A Model  of  Nonstaticnarity 

To  compute  the  properties  of  the  log  average  spectrum  and 
average  log  spectrum  for  nonstaticnary  processes  the  following  model 
is  proposed.  Ue  assume  that  the  nonstationary  process,  (x(t)I,  is 
produced  by  passing  ar  underlying  stationary  process,  iy(t)i, 
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through  a time-varying  linear  system.  Ue  further  assume  that  the 
linear  system  changes  slowly  enough  to  be  considered  stationary  over 
each  clata  segment  for  which  a periodogram  is  to  be  computed.  hus, 
the  resulting  process  mai  be  considered  stationary  for  each  segment. 
From  segment  to  segment,  however,  the  energies  at  each  frequency  may 
vary  considerably. 

The  statistical  fluctuations  in  the  process,  tnen,  are  modeled 
by  appropr iately  changing  the  frequency  response  of  this 
time-varying  linear  system.  Clearly,  these  changes  represent,  in 
some  sense,  a "schedule'"  of  nonstat  ionar  i ty.  In  general,  this 
schedule  is  a function  of  the  physical  phenomena  producing  the 
random  process.  One  would  expect,  for  example,  that  it  would  be 
approximately  the  same  for  different  recordings  of  the  same  musical 
selection  even  though  each  recording  may  represent  a statistically 
different  process. 

Admittedly  this  is  a simple  model.  Clearly,  many  signals  are 
nonstationary  over  very  short  intervals.  For  the  signals 
represented  in  figures  4.1  - 4.3,  the  individual  data  segments  were 
0.4096  seconds  long  whereas  speech  typically  exhibits 
nonstat ionari ty  over  much  shorter  intervals.  A model  encompassing 
such  detail,  however,  would  be  enormously  complex.  As  will  be  seen, 
a simpler  model  serves  well  in  understanding  the  log  spectral 
estimators. 

Ilathemat ical ly,  this  model  is  represented  by 

xk(t)  * b,(t)  ($  y,(t)  (4.1) 
where  x^tt)  is  a finite  sample  function  of  the  nonstationary  process 
x(t),  y,(t)  is  a finite  sample  function  of  the  underlying  stationary 
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process,  y(t),  and  b,(t)  is  the  impulse  response  of  the  ith  linear 
system.  The  relationship  in  (4.1)  is  only  approximately  equal  since 
the  convolution  is  performed  after  windowing  y(t).  In  practice,  it 
is  x(t)  that  is  windowed.  However,  by  taking  note  of  these  effects, 
we  may  proceed  with  the  analysis. 

Ue  now  compute  a periodogram  from  x,(t).  The  periodogram  is  a 
function  of  both  the  underlying  process  y(t)  and  the  ith  system 
response.  Thus,  we  have 

Ix„(f)  = (l/T)i:rix,(t))|J 

=*  (1/T)  IJ'iy.(t)©b,(t' I |J  ^ (l/T)  |Y,(f)-(3,(f)  |J 

* Iail(f)-iP,(f)iJ  (4.2) 

where  Y,(f),  £,(f)  are  the  finite  Fourier  transforms  of  y.(t)  and 

b,(t),  respectively,  and  Iw>,  ( f ) is  the  periodogram  associated  with 
Y(t).  Similarly,  the  loc  periodogram  is 

X x.t ( f ) =*  logIx.,(f)  2:  logIy„(f)  + logl^lf)!2  (4.3) 

and  the  smoothed  spectral  estimator  of  (2.55)  is 

p„(f)  = (i/N) z.^ix.jf)  = (i/fj) siN1iu.l(f)-ia,(f) i2  • (4.4) 

The  log  average  spectrum  and  average  log  spectrum  are  defined 
exactly  as  in  (3.2)  ano  (3.3),  but  are  written  in  terms  of  the 
periodogram  associated  with  y(t)  as 

Px(f)  - log  P»  ( f ) * log  ( (l/n)  2.?ilx..(f) ) 

2=  log[(l/M)  2iN1Ia„(f).|/3i(f)|2)  (4.5a) 

and 

L*  ( f ) - (un)  SiVoglxJf) 

* ( i/n)  2,% log Iy,, ( f)  + (l/N)  2.Nil°yl(Mf>l2 

* La(f)  + (i/n)  2,N1log1|3,if)l2  . (4.5b) 


Thus,  we  see  that  the  log  spectral  estimators  for  the  nonstationary 
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process  are  biased  from  an  estimate  of  the  log  spectrum  of  the 

underlying  stationary  process,  y(t),  by  amounts  that  are  a function 
of  the  time-varying  linear  system  and  the  method  of  averaging 

(linear  or  logari  thmic) . If  the  (3,(f)s  are  equal,  then  the  biases 
are  equal  and  represent  mererly  a gain  or  attenuation;  this,  of 
course,  is  expected  since  under  that  condition,  x ( t ) is  stationary. 

The  net  result  of  (4.5)  is  to  define  the  log  average  spectrum 

and  log  average  spectrum  for  x(t)  in  terms  of  the  log  spectral 

estimates  of  y(t).  The  proper  interpretation  of  y(t)  is  variable. 
The  precise  definition  of  y(t)  and  its  associated  log  spectrum,  as 
well  as  the  definition  of  (3,(f),  will  influence  the  biases 
associated  with  (4.5).  Nonetheless,  (4.5)  does  provide  a way  of 

comparing  and  understanding  the  log  spectral  estimators  of  x(t). 

4.3  Statistical  Properties 

In  computing  the  statistical  properties  of  P*(f)  and  Lx(f),  an 
important  point  should  be  kept  in  mind.  For  a given  experiment,  the 
|3,(f)s  are  deterministic.  For  fixed  i and  f,  Ia,,(f)  is  a random 
variable;  in  constrast,  |G,(f)  is  a constant.  Thus,  for  example, 

E !l(3,  ( f ) I1)  = 1/3,  ( f)  lJ  (4.6a) 

and 

var  (|(3, (f)  lJ!  = 0 . (4.6b) 

The  mean  and  variarce  of  Ix„(f)  and  P*(f)  are  easily  computed 
from  (4.2)  and  (4.4).  Proceeding  as  in  section  2.10  and  noting 
again  the  approximations  for  non-Gaussian,  non-white  processes,  we 
have 

E (I,„(f ) ) * l(3,(  f ) |**E  (!„.,  ( f ) 1 = 1(3,  ( f)  l:-Ga(f)  . (4.7a) 
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var  * 10. (f)  r-var  flaj  * 10.  ( f ) l4-Ga2<  f ) , (4.7b) 

E (P* ( f ) I a EilUN) 

= (i/N)2.NJ/3.(f)l2-EfIa.,(f)} 

fl!  Ga  ( f ) » (1/rj)  2.NJ0,(f>ls  (4.7c) 

and 

var(Px(f)!  a;  var  ( (i/rj)  ^I,,,  ( f ) ) 

= ((UN)sI.Nll/3l(f)l4)-var)Ii,.l(f)) 

= Gw* ( f ) - ( (l/NJ'I^l^tf)!4)  (4.7d) 

where  Ga(f)  is  the  spectral  density  function  of  the  stationary 
process  y(t).  Note  that  we  have  used  (4.G)  to  justify  bringing 
( f ) ( 2 out  of  the  expectation  and  variance  operators  as  a scaler. 

The  mean  and  variance  of  the  log  spectral  estimators  may  now  be 
derived.  This  is  easiest  for  the  average  log  spectrum  since  the 
effect  of  log  !& ( f ) l2s  is  additive.  From  (3.12)  and  (4.5),  and 
extending  (4.6)  to  log  1/3,  ( f)  I2,  we  have 

E IL*(  f ) ) a EIL,(f)J  + E((i/n)  I.^log  l(3,(f ) l*» 

= EIL„(f)l  + ( i/m)  S^Eiloglfl.mi’f 

a logGa(f)  - v + ( l/u)  log  1(3,  ( f ) I2  (4.8a) 

and 

var)Lx(f))  » var  !La  ( f ) 1 + var  { (i/n)  S^log  l(3,(f ) |M 

* n2/6N  2 1.8449—/N  . (4.8b) 

Corresponding  resuits  for  the  log  average  spectrum  are  somewhat 
more  difficult  and  require  an  additional  approximat ion.  In  general, 
the  (3,  ( f ) s are  different  for  each  i.  Thus,  Px(f)  is  now  a sum  of 
weighted  chi-square  randcm  variables.  Exact  computations  using  this 
distribution  yield  complicated  open  form  results  131) , (8,  p. 273)  . 
However,  as  was  done  in  section  2.9  for  S,(f),  we  approximate  P«(f) 
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as  a chi-square  random  variable  with  a constant  of  propor tional i ty 
and  degrees  of  freedom  given  by  (2.18).  Thus  Px(f)  = r*X2n  with 
n = EOF  (Px ( f ) ) = 2-EMPx(f)l  /var  f Px  ( f ) I 

- 2 (Gtf(f ) • I (l/N)  2^1/3.  (f)  !23  )2/  (G,J(f).[(UN)JI^I3i(f)l‘*]  ) 
= 2N-  C (l/N)  I^l^if)  I2]2/  ((UN)  (f)  I4]  (A.  9a) 

and 

r = E (P*  ( f ) ) / ri  = Ga  ( f ) • t ( lm)  2AI(3t(f)  I2]  /n  . (4.9b) 

Using  Table  2.1,  then, 

E(Px(f)l  - ^(n/2)  + log ( 2 r ) 

*=  ^(n/2)  - log(n/2)  + log  (Gu(f ) • ( (i/n)  XwJ&tf ) 1*1  ) 

* log Gy  ( f ) + log[(i/N)  2.Nil|3.(f)  I21  (4.10a) 

and 

var(Px(f))  * ^'(n/2)  * 2/n 

= C(l/N)  lAIP.tfJH  t (N-  [(l/N)  Z"iKMf  >!*]*>  • (4.10b) 

He  see,  then,  that  the  log  soectral  estimators  for  the 
nonstationary  process,  x(t),  are  biased  from  the  corresponding 
estimators  for  the  underlying  stationary  process  by  amounts  that  are 
a function  of  the  schedule  of  nonstationari ty . Since  this  schedule 
may  differ  from  process  to  process,  it  is  difficult  to  comment 
quantitatively  about  thei  effect  of  these  biases.  However,  some 
general  observations  are  possible. 

For  the  average  log  spectrum,  the  bias  is  the  logarithm  of  the 
geometric  mean  of  the  the  |(3,(f)l2s;  for  the  log  average  spectrum  it 
is  the  logarithm  of  the  arithmetic  mean.  As  discussed  in  section 
3.8,  the  geometric  mean  is  always  less  than  the  arithmetic  mean; 
thus  from  (3.9),  (un)  log  l(3t  ( f ) I2  < log  (l/N)  I0t  ( f ) I2.  The  way  in 
which  we  define  the  If3,(f)|:s  is  arbitrary  to  the  extent  thai  the 
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l(3v  ( f ) I’s  can  be  arbitrarily  scaled;  the  spectrum  associated  with 
the  underlying  process  is  then  adjusted  accordingly. 

Thus  the  biases  in  (4.10)  can  be  adjusted  to  be  completely 
negative  or  positive.  However,  the  ALS  bias  will  always  be 

numerically  less  than  the  LAS;  if  the  biases  are  negative,  the  ALS 
bias  will  be  more  negative  while  if  they  are  positve,  the  LAS  bias 
will  be  more  positive.  The  general  effect  is  that  peaks  will  tend 
to  be  accentuated  in  the  log  average  spectrum  while  troughs  will  be 
accentuated  in  the  averye  iog  spectrum. 

If  the  0,(f)s  are  the  same,  the  bias  terms  are  equal  and  act 
merely  as  scalers.  If  they  specifically  equal  one,  then  (4.S)  and 
(4.10)  reduce  to  the  results  for  stationary  processes,  (3.11)  and 
(3.12) . 

Both  estimators  are  consistent.  Interestingly,  while  the  ALS 
variance  is  the  same  as  for  a stationary  signal  the  LAS  variance  is 
a function  of  the  nonstationari ty.  For  particular  values  of  the 
(3,  (f)s,  it  may  even  be  greater  than  the  ALS  variance.  In  fact,  this 
i3  what  we  observe  in  figures  4.1  - 4.5. 

(lore  quantitative  insight  into  (4.8)  and  (4.10)  is  possible  by 
extending  our  model  somewhat.  Specifically,  note  that  in  these 

equations,  (l/n)  2,NJ0,  ( f ) Is  is  a sample  mean  of  10(f)  I*  (computed  from 
the  N samples,  0,(f));  (i/n)  log  10.  < * > I*  is  a sample  mean  of 

log  10  ( f ) I’;  and  so  on.  By  assuming  a particular  form  for  the 
distribution  of  the  10,(0  1’  from  segment  to  segment,  we  can  make 
more  precise  comments  about  the  relationship  of  the  log  spectral 
estimators. 

In  light  of  the  above  comments,  we  first  rewrite  (3.8)  and 
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(3.10)  in  terms  of  the  aopropriate  expectations  of  l(3(f)l*.  If  N is 
large  enough,  this  is  a reasonable  approximation.  Thus  we  have 

E (L*(f ) ) * leg Gy ( f ) — y + E 1 log  1(2 ( f ) I21  . (4.11a) 

var(Lx(f)>  a n!/6N  = 1.6449... /N,  (4.11b) 

EfMfH  a logGy  ( f ) + logE  { 1/3  ( f ) I2)  , (4.11c) 


and 

var  {py(  f ) I ^ E (l(3(f)  I4!  / (N-E2  (|(3(f)  I2)  ) . (4. lid) 

There  are  several  reasonable  choices  for  the  distribution  of 
l(3(f)|J.  Undoubtedly  the  distribution  is  different  for  each 

experiment.  For  the  purposes  of  this  research,  however,  we  assume 
that  the  logarithm  of  l|3(f)|2  is  uniformly  distributed  on  some 
interval,  (a.b).  This  is  equivalent  to  stating  that  the  logarithm 
of  the  power  in  the  time  varying  system  is  linearly  distributed. 
Although  possibly  a cruce  model,  it  is  not  unreasonable  for  vocal 
and  instrumental  signals  since  it  is  well  known  that  the 
conventional  musical  nolation;  pp,  p,  mp,  mf,  f,  ff;  represents 
logarithmically  increasing  amplitudes.  Also,  experimental  results 
(described  in  section  4.6)  suggest  that  it  is  close  enough  to  allow 
meaningful  results  to  be  derived  that  are  helpful  in  gaining  a more 
intuitive  understanding  of  log  spectral  estimation  for  nonstationary 
processes. 

Uith  this  assumption,  the  results  of  Appendix  B,  section  B.3, 
are  applicable.  If  we  let  X = logj/3(f)l2  and  Y = 1/3  ( f ) 1 2 then 

(B.29)  and  (B.30)  can  be  substituted  in  (4.11).  Note  that  a < 


1 og  10(f)  I*  < b,  A < I /3 ( f . I*  < B.  0 « B/A  is  the  dynamic  range  of 
|(3(f)l2,  and  d = log(D).  Note  also  that  the  dynamic  range  of  l|2(f)l2 
will  in  general  c • a function  of  frequency,  f.  However,  this  is  not 
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shown  explicitly.  Rewriting  (4.11),  then,  we  have 

E 1LX( f ) ) * logGy (f ) - Y + b - d/2,  (4.12a) 

var  1L* ( f ) ) » n'/SN  » 1.G449—/N,  (4.12b) 

E ( pK ( f ) 1 as  logGa(f)  + b - log(d)  + log(l  - 1/D) 

a:  logGa(f)  + log(B)  - log(d),  (4.12c) 


and 

var  ( pK ( f ) 1 * (3*/2d)  / CN»  (BI/dJ)  ] *d/2N  . (4.12d) 

From  (4.12)  we  can  see  that  for  this  model,  the  bias  terms  are 
a function  of  the  dynamic  range  of  13(f)  I2.  The  +b  terms  merely 
scale  the  estimators  by  the  same  amount.  Since,  as  d increases,  d/2 
increases  much  more  rapidly  than  log(d),  the  ALS  estimator  will 
generally  be  influenced  more  by  the  nonstationarity. 

Moreover , 

var  1L*( f ) ) / var-  ip* ( f ) 1 a:  it2 / 3d  = 3. 2938” • / d (4.13) 

and  we  see  that  as  the  dynamic  range  gets  larger , this  ratio  of 
variance  gets  smaller.  If  d > 3.299S’"  the  log  average  spectrum 
will  have  a larger  variance  and  will  thus  be  less  stable.  This 
value  corresponds  to  d > 14.3  dB.  As  shown  in  section  4.6,  for  the 
data  of  figure  4.1  (with  logl3(f)l2  uniformly  distributed)  d ^ 
IS. 89  & 69.0  dB.  Hence,  one  would  expect  the  LAS  variance  to  be 
about  4 times  that  of  the  ALS,  or,  equivalently,  the  LAS  standard 
deviation  to  be  twice  that  of  the  ALS.  Although  difficult  to  judge, 
this  seems  to  be  in  reasonable  agreement  with  figure  4.1. 

During  the  course  of  this  research,  an  interesting  collateral 
issue  arose.  The  results  expressed  in  (4.10)  are  derived  by 
assuming  the  distribution  of  P„(f)  to  be  proportional  to  a 
chi-square  distribution  with  the  degrees  of  freedom  given  by  (2.18). 
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As  an  alternative  approach,  the  same  calculations  were  made  by 
approximating  Px  { f ) with  a lognormal  distribution  (2.3).  The 
motivation  for  this  is  the  observation  that  logPx(f)  has  a 
distribution  which  tends  to  be  more  Gaussian  than  Px ( f ) ; this  is 
particularly  true  for  large  N.  Thus,  by  assuming  Px(f)  to'  be 
lognormal  and  px(f)  to  be  normal  the  desired  statistics  of  px(f)  can 
be  derived. 

Using  this  approach,  the  mean  and  variance  of  Px(f) 

corresponding  to  (4.10)  ere 

E(px(f)l  * 1 og  Gu  ( f ) + log  [ (i/n)  ( f ) H 

- ( 1/2 ) var  { px ( f ) 1 (4.14a) 

and 


var(px(f)l  * logtl  + 

+ Him)  I.?J(3,(f)r)  / (N[(i/N)ZiMitfMf)l,l*H 
* ((1/w)  INil(3.(f)l*)  / (N  £ (l/N)  Xi?J0Jf)l*]*)  (4.14b) 

It  is  interesting  that  (4.10)  and  (4.14)  are  so  similar.  The 

approximating  expressior s for  the  variance  are  the  same;  the 
expressions  for  the  mean  are  the  same  if  the  digamma  function  in 
(4.10)  is  approximated  by  the  first  two  terms  of  its  series 
expansion  (A. 11 ) . 

Computations  using  (4.14)  agree  well  with  the  data  of  section 
4.6.  For  small  N,  the  lognormal  approach  tends  to  predict  values 
higher  than  the  observed  mean  and  lower  than  the  observed  variance. 
The  chi-square  approach  coes  just  the  opposite;  however,  the  results 
generally  agree  more  closely.  Since  the  results  of  this  alternative 
approach  yield  no  new  information  (beyond  the  fact  that  two 
seemingly  unrelated  distributions  give  such  similar  results). 
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details  of  its  derivatior  are  not  given. 

4.4  The  Activation  Spectrum 

He  again  define  the  activation  spectrum,  A*(f),  as  the 
difference  between  the  average  log  spectrum  and  the  log  average 
spectrum  (3.7).  The  expected  value  of  A*(f)  is  derived  from  (4.8a) 
and  (4.10a)  as 

E (Ax  ( f ) 1 as  -y 

+ ( (i/N)  I ^loglfMf)!*  - log[(i/N)  Z^l(3.(f)l2) ) . (4.15) 

In  terms  of  the  assumed  uniform  distribution  of  log|$(f)|J,  this 
becomes 

E(A„(f)}  a -y  +•  (E  t log  1(3  i f ) I2)  - logE  (l(3(f ) 1*1  ) 

a;  -Y  - Cc/2  - log(d))  . (4.1G) 

From  (4.15)  we  see  that  the  activation  spectrum  has  a value 
equal  to  or  less  than  -y.  From  (4.15)  we  see  that  this  additional 
amount  is  the  logarithm  cf  the  ratio  of  the  geometric  and  arithmetic 
means  of  the  l(3,(f)lJs.  Thus  the  lower  bound  of  this  ratio, 
developed  by  Cox  (3.10),  is  applicable.  As  the  dynamic  range  of  the 
nonstationarity  decreases,  the  term  in  parenthesis  in  (4.15)  becomes 
less  negative  and  (4.15)  approaches  (3.13),,  the  result  for  a 
stationary  process.  The  net  result  is  that  the  greater  the 
nonstationary,  represented  by  an  increased  dynamic  range  of  the 
l(3,(f)lls,  the  more  negative  Ax(f)  will  tend  to  be. 

This  is  explicitly  clear  in  (4.1G).  In  this  equation,  the  term 
in  brackets  is  given  as.  a direct  function  of  dynamic  range.  Since 
this  function  is  monotonical ly  decreasing  for  increasing  d (as 
easily  seen  by  computing  its  derivitive),  the  greater  the  dynamic 
range  of  the  l/3,(f))Js  the  more  negative  the  activation  spectrum. 
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This  result  is  perhaps  one  of  the  most  useful  of  this  research. 
The  activation  spectrum  provides  a sensitive  test  for 
nonstationarity  and  gives  a distribution  of  nonstationari ty  with 
frequency.  This  is  the  motivation  for  terming  Ax(f)  the  activation 
spectrum;  practical  signals  often  consist  of  a nonstationary  process 
plus  stationary  noise.  An  examination  of  Ax(f)  for  such  a process 
reveals  those  frequencies  in  which  the  activity  of  the  signal 
predominates. 


Turning 

again 

to  figures  4.1  - 4.5, 

inspection 

of  the 

activation 

spectra 

provides  insight  into 

the  nature 

of  the 

represented 

signals. 

In  figure  4.1,  for 

example,  Ax(f) 

has  a 

characteristic  "necKlace’  shape.  This  results  from  (1)  the  fact 
that  this  recording  of  Caruso  has  a large,  resonant  peak  in  the 
spectrum  near  700  Hz  and  (2)  the  strong  presence  of  stationary 
surface  noise.  In  the  vicinity  of  this  peak,  the  signal 

predominates  and  Ax(f)  has  a deep  trough.  On  either  side  of  700  Hz, 
however,  the  S/N  ratio  slowly  decreases.  The  result  is  a gradual 

lessening  of  the  nonstat. ionary  character  of  the  signal  indicated  by 
a positive  increase  in  Ax(f).  Outside  the  effective  bandwidth  of 
the  singing  (160  Hz,  3250  Hz)  Ay(f)  is  close  to  -2.5  dB  indicating 
that  the  signal  is  essentially  stationary  noise.  It  is  interesting 
to  note  that  listening  tc  this  recording  through  a sharp  cut-off  low 
pass  filter  indicates  that  there  is  no  audible  music  energy  below 
160  Hz,  the  cutoff  indicated  by  Ax(f). 

Similar  observations  may  be  made  about  the  other  figures.  The 
activation  spectrum  of  the  female  singer  (figure  4.4)  is  very 
negative  indicating  the  extensive  dynamics  of  the  singing.  It  is 
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also  apparent  that  there  is  little  activity  in  the  lower 
frequencies.  Conversely,  the  activation  spectrum  of  the  string 
ensemble  (figure  4.5)  shows  the  presence  of  activity  at  both  very 
low  and  very  high  frequencies.  This,  of  course,  would  be  expected 
from  the  uide-band  nature  of  musical  instrumentation. 

Another  interesting  point  is  the  apparent  correlation  between 
the  spectral  dynamic  ranee  of  a process  and  the  activation  spectrum. 
Again,  this  is  evident  in  the  figures.  Clearly  the  spectral  dynamic 
range  is  influenced  not  only  by  the  randomness  of  the  periodograms. 


but  by 

the  1(3.  ( f ) l2s. 

As  their 

dynamic 

range 

increases , 

the 

spectral 

dynamic  range 

increases 

and,  as 

shown 

by  (4. Id) , 

the 

activation  spectrum  beccmes  more  negative.  This  is  particularly 
interesting  since  the  work  of  Cox  predicts  only  that  the  activation 
spectrum  is  bounded  by  a curve  that  is  a function  of  the  spectral 
dynamic  range.  These  results  show  that  not  only  is  Ax(f)  bounded  by 
this  curve,  but  will  generally  be  statistically  correlated  with  it 
(and,  thus,  to  the  spectral  dynamic  range). 

One  other  phenomenon  present  in  these  figures  is  the  presence 
of  occasional  peaks  in  the  activation  spectrum  with  values  greater 
than  -2.5  dB  (occasionally  very  close  to  zero).  These  represent  the 
presence  of  a coherent  component  of  the  signal.  Since  such  a 
component  does  not  statistically  vary,  Ax ( f ) is  zero  (although  the 
presence  of  some  noise  ir  these  components  prevents  it  from  actually 
being  zero).  In  these  particular  figures,  the  peaks  are  most 
generally  at  60  Hz  or  120  Hz  representing  sinusoidal  hum  in  the 
electronics  used  for  reproducing  and  recording  the  signals. 
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4.5  The  Effect  of  Additive,  Stationary  Noise 

Up  to  this  point,  our  derivations  and  discussions  have  been  in 
terms  of  noiseless  signals.  In  the  case  of  stationary  processes, 
this  is  not  significant  since  noise  is  (usually)  stationary  and 
additive;  the  resulting  process  is  therefore  also  stationary.  Noise 
is  of  consequence  only  if  it  is  desired  to  separate  its  spectrum 
from  that  of  the  signal. 

For  nonstationary  processes,  however,  the  effect  can  be 
significant.  Uhere  the  noise  is  large  compared  to  the  signal,  the 
resulting  process  behaves  as  a stationary  process;  where  the  noise 
is  comparati vely  small,  the  signal  acts  as  a nonstationary  process. 
Thus  it  is  important  to  ciscuss  its  effect  further. 

It  is  easiest  to  understand  the  effect  of  additive  noise+  if  we 
assume  that  it  alters  the  distribution  of  the  1/3^  ( f ) l2s  and,  thus, 
is  absorbed  in  (4.15).  Other  approaches  would  be  considerably  more 
di f f icul t . 

To  understand  just  how  noise  perturbs  the  lf3,(f)|Js,  consider, 
first,  noise  added  to  a stationary  process.  If  the  signal  and  noise 
are  uncorrelated,  as  is  usually  the  case,  and  letting  x(t)  *=  y C t ) + 
n(t),  then  G„ ( f ) - Gy(t)  + G„(t)  where  x(t),  y(t)  are  stationary 
processes  and  n(t)  is  stationary  noise.  Then  logG*(f)  = 

1 og ( Gu ( f ) + G„(f)).  Clearly,  log  G„ ( f ) is  bounded  below  by  the 

larger  of  Gv ( f ) and  GN(f). 

In  the  case  of  nonstationarity,  as  modeled  in  section  4.2,  the 


+Not  all  noise  in  a signal,  of  course,  is  necessarily  additive. 
Many  situations,  however,  can  be  accurately  modeled  as  such. 
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effect  of  the  nonstation. arity  is  to  reduce  or  increase  the  spectrum 
of  the  underlying  stationary  process,  y(t),  by  an  amount  dependent 


upon  each  l(3,(f)|J. 

If  noise, 

however, 

has 

been 

added  to 

the 

nonstationary  process 

, then  no 

matter 

how 

srnal  1 

a particular 

l(3,(f)!J,  the  spectrum 

will  not 

be  less 

than 

the 

value  of 

the 

spectrum  of  the  noise.  The  effective  dynamic  range  of  the  ( f ) I *s 
is  thus  reduced.  The  noise  has  effectively  raised  their  minimum 
possible  value. 

As  the  noise  increases,  the  interval  t A , B 3 over  which  l (3  ( f ) 1 2 
is  distributed  decreases  while  A increases.  Simultaneously, 
(i/n)  log  IP.  ( f > I2  ancJ  log  [ ( i/n)  ( f ) |J3  increase  until,  when  the 

noise  is  great  enough,  they  are  equal,  their  difference  is  zero  and 
the  bias  terms  in  (4.8)  and  (4.10)  now  represent  the  spectrum  of  the 
noise.  Accordingly,  E I A,  ( f ) } will  now  be  -2.5  dB  as  expected  for  a 
stationary  process. 

Since  u/N)  £.Ni  log  1(3,  ( f ) I*  < log  (i/n)  13.  ( f ) I*  and  both  quantities 

increase  -with  the  addition  of  noise,  the  first  quantity,  the  ALS 
bias,  will  be  affected  the  most.  Thus,  the  average  log  spectrum  is 
most  sensitive  to  the  adcition  of  noise.  One  can,  in  fact,  conceive 
of  situations  where,  i~  the  difference  in  the  ALS  and  LAS  is 
significant  enough  and  the  noise  large  enough,  the  average  log 
spectrum  would  be  almost  completely  engulfed  by  the  spectrum  of  the 
noise.  This  is  an  important  consideration  in  deciding  which 
estimator  to  use  in  a particular  situation.  As  will  be  seen  in 
chapter  5,  this  has  important  consequences  when  using  log  spectral 
estimators  to  estimate  linear  system  functions. 
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4.6  Experimental  Results 

Figure  4.6  represents  the  result  of  computing  log  spectral 
estimates  for  a simulated  nonstationary,  Gaussian,  white  process. 
This  process  was  generated  by  multiplying  each  data  segment  by  a 
random  number  with  a hyperbolic  distribution  (see  Appendix  C) . This 
is  equivalent  to  the  ime-varying  linear  system  discussed  in  section 
4.2  being  a time-varying  amplifier.  For  reasons  to  be  discussed 
later  in  this  section,  the  distribution  of  l$(f)  lJ  was  selected  to 
have  a dynamic  range  of  d = 15.388  = 69.0  dB  over  the  inerval  LA.B1 
where  A = 0.00075377  and  B - 6035.05.  For  these  values,  ElLx(f))  ^ 
30.8  dB,  E (px(f ) ) * 55.3  dB,  and  E)AK(f))  & -25.0  dB.  Also,  (4.13) 
predicts  that  var(Px(f)l  z 4.8-var  )px(f)  1 (about  twice  the  standard 
deviation).  All  these  values  are  in  apparent  agreement  with  figure 
4.6. 

Tables  4.1  - 4.4  are  the  sample  means  and  sample  variances  for 
the  data  of  figure  4.6  compared  to  the  theoretically  predicted 
values.  These  tables  are  similar  to  tables  3.1  - 3.4.  As  before, 
agreement  between  theory  and  the  empirical  result  is  good. 

Figure  4.7  presents  the  results  of  a more  significant 
simulation  experiment.  In  this  case,  the  nonstationary,  white 
process  of  figure  4,6  was  passed  through  two  linear  systems  to  color 
the  process,  and  then  added  to  stationary  noise.  The  particular 
system  and  amount  of  noise  was  chosen  in  an  attempt  to 
parametrically  duplicate  the  real  process  depicted  in  figure  4.1 
(Caru90  9inging  from  an  cld  acoustic  recording). 

From  figure  4.1,  we  note  that  the  maximum  displacement  in  the 
activation  spectrum  is  approximately  -25.0  dB.  Utilizing  our  model 


FIGURE  4.6 

Log  spectral  estimates  (N  470)  for  a simulated  nonstationary, 
Gaussian,  white  process  using  models  proposed  in  the  text:  (a)  log 
average  spectrum,  (b)  average  log  spectrum,  (c)  activation  spectrum, 
and  (d)  spectral  dynamic  range. 
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TABLE  4.1 

SAMPLE  MEAN  CF  THE  LOG  AVERAGE  SPECTRUM  FOR  A 
NONSTATIONARY,  GAUSSIAN,  WHITE  PROCESS  ~ N(O.IOOO) 


N 

Lower 

Sample 

Upper 

Predict 

Conf id. 

Mean 

Conf id. 

Value 

1 

55.3390 

57.08G9 

58.8348 

57.0616 

2 

53.9090 

54.0772 

54.2454 

54.1206 

4 

52.0933 

52.2284 

52.3585 

51.9698 

6 

5G.8090 

5G.9198 

57.0306 

56.8528 

8 

55.G504 

55.7587 

55.8670 

55.6792 

10 

54.G814 

54,7897 

54.8980 

54.7102 

20 

55.9079 

55.9733 

56.0387 

55.9442 

30 

57.2617 

57.3150 

57.3683 

57.3414 

40 

5G.5273 

56.5759 

56.6245 

56.6199 

50 

55.5GG4 

55.6143 

55.6634 

55.6587 

100 

55.57G2 

55.6115 

55.6468 

55.6360 

200 

5S.G499 

55.6767 

55.7035 

55.6977 

400 

5G.0007 

56.0203 

56.0399 

56.0203 

TABLE  4.2 


SAMPLE  MEAN  CF  THE  AVERAGE  LOG  SPECTRUM  FOR  A 
NONSTATIONARY,  GAUSSIAN,  WHITE  PROCESS  ~ N (0,1000) 


N 

Lower 

Sample 

Upper 

Predict 

Conf id. 

Mean 

Conf id. 

Value 

1 

55.3390 

57.0863 

58.8348 

57.0616 

2 

32.2972 

32.4165 

32.5358 

32.3380 

4 

33.4278 

33.5127 

33.5976 

33.4560 

6 

39.4829 

39.5545 

39.6251 

39.5445 

8 

37.5603 

37.6209 

37.6815 

37.5776 

10 

31.7308 

31.7860 

31.8412 

31.7455 

20 

34.3654 

34 . 4044 

34.4434 

34.3593 

30 

37.4663 

37.4974 

37.5285 

37.4699 

40 

35.3651 

35.3924 

35.4197 

35.3715 

50 

32.0416 

32.0661 

32.0906 

32.0146 

100 

31.7146 

31.7325 

31.7504 

31.6973 

200 

30.6875 

30.7003 

30.7141 

30.6549 

400 

30.4686 

30.4788 

30.4890 

30.4321 

7G 


TABLE  4.3 


SAMPLE  VARIANCE  OF  THE  LOG  AVERAGE  SPECTRUM  FOR  A 
NONSTATIONARY.  GAUSSIAN,  WHITE  PROCESS  ~ N (O.IOOO) 


N 

Lower 
Conf id. 

Sample 

Mean 

Upper 
Conf id. 

Predicti 

Value 

1 

28.9323 

30.1995 

31.5515 

31.0254 

2 

28.9058 

30.1712 

31.5220 

27.2819 

4 

17.278G 

18.0350 

18.8424 

23.3037 

6 

12.7498 

13.0793 

13.9038 

15.6508 

8 

11.9725 

12.49SG 

13.0561 

15.2449 

10 

11.9724 

12.4965 

13.0560 

15.2446 

20 

4.3713 

4».  5627 

4.7669 

4.8791 

30 

2.9015 

3.0286 

3.1642 

3.4021 

40 

2.4137 

2.5194 

2.6322 

2.8277 

SO 

2.4042 

2.5094 

2.6218 

2.8182 

100 

1.2705 

1.3261 

1.3855 

1.4421 

20» 

0.7349 

0.7671 

0.8014 

0.7437 

400 

0.3914 

0.4085 

0.4268 

0.3926 

TABLE  4.4 

SAMPLE  VARIANCE  OF  THE  AVERAGE 

LOG  SPECTRUM  FOR  A 

NONSTATIONARY, 

GAUSSIAN,  WHITE 

PROCESS  -v 

N (0, 1000) 

N 

Lower 

Sample 

Upper 

Predic  ti 

Con  fid. 

Mean 

Conf id. 

Value 

1 

28.9329 

30.1995 

31.5515 

31.0254 

2 

14.SA32 

15.1798 

15.8594 

15.5127 

4 

7.3668 

7.6893 

8.0335 

7.7564 

6 

5.0955 

5.3186 

5.5567 

5.1709 

8 

3.7566 

3.9210 

4.0966 

3.8782 

10 

3.1157 

3.2521 

3.3977 

3. 1025 

20 

1.5516 

1.6196 

1.6921 

1.5513 

30 

0.9867 

1.0299 

1.0760 

1.0342 

40 

0.7586 

0.7913 

0.8272 

0.7756 

50 

0.6149 

0.6418 

0.6705 

0.6205 

100 

0.3265 

0.3408 

0.3561 

0.3103 

200 

0.1818 

0.1898 

0.1983 

0.1551 

400 

0.1059 

0.1105 

0.1155 

0.0776 
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of  10(f)  lJ,  this  corresponds  to  a dynamic  range  of  d =*  15.888  a 59.0 
dB.  This  value  is  obtaired  by  numerically  solving  (4.15)  for  d.  By 
noting  the  value  of  the  log  average  spectrum  in  figure  4.1(a),  and 
utilizing  the  fact  that  the  underlying  stationary  process  is 
distributed  as  N(0,1000),  the  values  of  A and  B given  earlier  were 
also  derived.  These  parameters  were  then  used,  as  described  in 
Appendix  C,  to  produce  the  process  depicted  in  figure  4.6 

This  nonstationary  process  was  then  colored  by  passing  it 
through  two  linear  systems.  The  frequency  response  of  the  first, 
figure  4.8(a),  was  derived  from  a linear  average  of  the  log  average 
and  average  log  spectra  of  figure  4.2.  This  system  has  the  effect 
of  coloring  the  spectrum  of  the  simulated  process  to  approximate 
that  of  a typical  spectrum  of  singing.  In  producing  this  system, 
the  log  spectral  estimates  were  further  smoothed  by  convolving  them 
with  a spectral  window.  Note  that  this  is  similar  to  the  frequency 
smoothing  discussed  in  section  2.9;  in  this  case,  however,  it  is  the 
log  spectral  estimates  that  are  smoothed. 

The  process  was  ther  filtered  again  with  the  system  depicted  in 
figure  4.8(b).  This  is  a linear  combination  of  the  LAS  and  ALS 
estimates  (see  chapter  5)  of  the  frequency  response  of  the  acoustic 
recording  horn  that  produced  the  sharp  resonant  peaks  in  the  data  of 
figure  4.1.  This  was  to  simulate  those  sharp  resonances. 

The  additive  noise  was  generated  by  passing  a stationary  random 
process  (see  figure  3.1)  through  the  system  shown  in  figure  4.8(c). 
The  shape  of  this  freqency  response  was  produced  from  log  spectral 
estimates  computed  from  passages  of  the  Caruso  recording  containing 
surface  noise  only;  as  such,  it  is  an  estimte  of  the  spectrum  of 
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that  noise.  The  value  of  this  simulated  "Caruso"  noise  was  then 
scaled  to  a value  representative  of  actual  surface  noise  and  added 
to  the  simulated  nonstal ionary  process.  To  add  a final  touch  of 
comparability,  a GO  Hz  sinusoid  was  added  to  simulate  a coherent 
peak. 

As  can  be  seen  from  figure  A. 7,  the  resulting  log  spectral 
estimates  are  quite  comparable  to  figure  4.1.  Note  particularly 
that  the  activation  spectrum  exhibits  the  same  charac ter i s t ic 
"necklace"  shape.  This  correlation  demonstrates  that  the  models 
proposed  in  this  chapter  do  a reasonable  job  of  explaining  the 
marked  differences  observed  in  log  spectral  estimates  of 
nonstationary  processes. 
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FIGURE  4.7 

Log  spectral  estimates  (N  = 470)  for  a simulated  nonstationary, 
Gaussian,  colored  process  with  additive,  stationary,  colored  noise. 
Compare  to  figure  4.1:  (a)  log  average  spectrum,  (b)  average  log  spec- 
trum, (c)  activation  sDectrum,  and  (d)  spectral  dynamic  range. 
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FIGURE  4.8 

System  frequency  responses  used  in  production  of  the  simulated 
nonstationary  process  of  figure  4.7:  (a)  modern  recording  of  singing 
from  fiaure  4.2,  (b)  resonant  response  of  the  recording  horn  from  the 
data  of  figure  4.1  (see  chapter  5),  (c)  stationary  surface  noise  from 
the  data  of  figure  4.1. 
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CHAPTER  5 

AN  APPLICATION  OF  LOG  SPECTRAL  ESTIMATORS 

5.1  Digital  Log  Spectral  Estimation 

Llith  the  advent  of  modern  computer  technology,  computation  of 
spectral  estimates  has  become  a practical  reality.  High  speed 
techniques  enable  rapid  computation  of  Fourier  transforms  and,  thus, 
periodograms.  Similarly,  high-speed  convolution  [321  and  digital 
filter  design  enable  practical  spectral  smoothing  and  its 
application  to  linear  systems. 

As  mentioned  previously,  a direct  app.ication  of  log  spectra  is 
to  the  estimation  of  linear  system  functions.  In  the  remainder  of 
this  chapter,  we  will  discuss  this  application  in  the  context  of 
both  log  average  and  average  log  spectra. 

5.2  Blind  Deconvolution 

In  practice,  signals  are  frequently  encountered  that  are  the 
convolution  of  two  other  signals.  For  example,  a blurred  photograph 


is  the  convolution 

of  an  image 

with  a 

point 

spread 

function 

representing  the  out-of-focus 

or 

moving 

lens. 

Other 

similar 

situations  arise  in 

acoustics, 

geophysics, 

etc. 

The  problem  of 

separating  such  signals  is  called  deconvolution  and  is  the  topic  of 
much  current  research  (e.g.,  see  [331,  [341,  and  [351). 

In  some  instances,  one  of  the  signals  is  Known  and  it  is  a 


straight  forward  matter  to  recover  the  other. 


However,  a more 
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complex  problem  is  to  separate  them  when  both  are  unknown;  this 
problem  has  come  to  be  known  as  blind  deconvolution  [15]  . t"  This 
I more  difficult  problem  is  simplified  if  one  of  the  unknown  signals 

ie  of  a smaller  extent  than  the  other  (as  is  often  the  case  when  a 
long  selection  of  speech  or  singing  is  passed  through  a linear, 
stationary  system).  In  this  situation,  the  different  extents 
provide  a distinguishing  characteristic  needed  to  separate  the 
signals. 

Ue  will  find  it  convenient  to  proceed  in  terms  of  the  specific 
situation  found  in  old,  acoustic  recordings.  These  were  often  made 
using  the  recording  apparatus  depicted  in  figure  5.1  [151.  The 
musical  signal  was  amplified  by  an  acoustic  horn  which  in  turn  drove 
a stylus  as  it  cut  a grove  in  a wax  disc.  As  indicated,  the  singing 
signal,  s(t),  appearing  at  the  mouth  of  the  horn  was  affected  by 
passage  through  the  horn.  Speci f ical  ly,  it  was  convolved  with  the 
impulse  response  of  the  horn  and  the  resulting  recorded  signal, 
v(t),  was  badly  resonated.  Further  degredation  in  the  playback 
signal,  p(t),  resulted  from  addtive  surface  noise. 

* If  h(t)  is  known,  a restoration  filter  can  be  derived  and  the 

f- 

»E‘; 

signal  deresonated.  However,  h(t)  generally  varied  from  recording 
to  recording  and,  thus,  must  be  estimated  a priori.  Ule  will  first 
derive  a solution  to  this  problem  in  terms  of  the  homomorphic  theory 
of  Oppenheim  [16]  and  ils  relationship  to  log  spectral  estimation 


+ The  discussion  of  blind  deconvolution  and  its  specific 
application  to  deresoncting  acoustic  recordings  is  discussed  more 
fully  in  C15]  which  was  co-authored  by  the  author  of  this  document. 
T he  reader  is  refer ed  to  that  article  to  comp liment  the  brief 
discussion  presented  here. 
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via  the  average  log  spectrum.  This  leads  naturally  to  a similar 
approach  using  the  log  average  spectrum. 

5.3  Homomorphic  Deconvolution 

The  theory  of  homomorphic  filtering,  or  generalized  linear 
filtering,  is  an  extension  of  the  familiar  linear  system  theory.  A 
linear  system  has  the  characte- istic  that  it  obeys  superposition 
across  addition.  If  the  response  of  the  system  to  x,(t)  is  y,(t). 
and  the  response  to  x?(t)  is  y..(t),  then  the  response  to  r,x4(t)  + 
r?x?  ( t ) is  r1y1(t)  + r?y,(t).  A generalized  linear  (or  homomorphic) 
system  has  the  same  property  but  across  one  of  a broad  class  of 
operators. 

Sucfj  a system  transforms  the  signals  so  that  the  desired 
operation  is  mapped  into  addition.  At  this  stage,  linear  filtering 
may  be  introduced  into  the  system.  For  example,  if  two  signals  have 
been  convolved,  a Fourier  transform  (mapping  convolution  into 
multiplication)  followed  by  a complex  logarithm  (mapping 

mu  1 1 ip  1 ica ton  into  addition)  results  in  the  transformed  signal  being 
the  sum  of  the  two  inputs.  After  linear  filtering,  the  inverse 
logarithm  and  Fourier  transform  are  performed.  The  result  is  that 
if  x,(t)®x,(t)  is  an  input  to  the  system,  then  the  output  is 
y,  ( t)  ®y?(  t) . 

The  homomorphic  theory  can  be  naturally  applied  to  blind 
deconvolution.  For  example,  the  frequency  response  of  a system 

through  which  a signal  had  been  processed  might  be  estimable  if  the 
resulting  convolved  signal  were  mapped  into  a sum  of  the  component 
signals  by  a homomorphic  syster  as  described  above.  This  estimate 
could  then  be  used  to  prcduce  a restoration  filter. 
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Proceeding  according  to  this  theory,  then,  we  express  the 
relationships  of  figure  5.1  mathematically  as 


p(t)  = 

s ( t ) ® h ( t ) + 

n ( t) 

(5.1) 

where 

p(t)  is  the  playback 

signal,  s(t)  the  singing, 

h(t) 

the 

resonant  impulse 

response  and 

n(t)  the  additive 

noise. 

At 

this 

point. 

we  assume 

the  process  to  be  noiseless 

(we  consider 

the 

effects  of  noise  later).  Taking  the  Fourier  transform  and  complex 
logarithm  of  both  sides  of  (5.1),  we  have 

logP(f)  = logS(f)  + logH(f)  (5.2) 

At  this  stage,  one  might  consider  averaging  several  recordings 

(members  of  the  ensemble)  so  that  the  effect  of  logS(f)  would 

reduce  to  zero  (or  some  constant)  leaving  an  estimate  of  logH(f). 
Note  here  that  h(t)  is  deterministic  while  we  consider  s(t)  to  be 
random.  Unfortunately,  only  one  recording  is  usually  available. 
Therefore,  following  the  procedure  of  earlier  chapters,  we  assume 
the  process  to  be  ergodic,  and  average  over  adjacent  data  segments. 
Thus  we  have 

p,(t)  = wt  ( t ) • p ( t ) = w,(  t)  • ts  ( t)  ®h  ( t)  ] 

* s,(t)®h(t)  . (5.3) 

The  approximate  equality  in  (5.3)  results  from  a consideration  of 
windowing  effects.  If  the  windows  are  long  and  smooth  enough,  then 
the  approximation  holds  closely. 

Applying  the  Fourier  transform  and  complex  logarithm  once 
again,  (5.3)  becomes 

logPt(f)  * logS,(f)  + logH(f)  . (5. A) 

Since  this  a complex  logarithm,  we  can  rewrite  (5. A)  in  terms  o*  •• 
real  and  imaginary  parts.  Doing  this  and  averaging  over  the  N 


8G 


segments  gives 

(l/N)  ^i^log  IPJf)  | * 1 og  IH ( f ) I + (i/n)  ^T^log  ISk(f ) I (5.5a) 

and 

(1/N)  ZiVPilf)  * £H(f)  + (1/N)  . (5.5b) 

Ue  desire  that  the  last  term  in  (5.5a)  and  (5.5b)  will  average 
to  zero  leaving  an  estimate  of  the  frequency  response  of  the 
resonating  system.  In  general,  however,  this  does  not  happen.  This 
is  because  the  last  term  in  (5.5a)  is  really  in  the  form  of  a log 
average  spectrum  and,  as  such,  is  an  estimate  of  the  spectrum  of 
s(t).  Clearly,  such  a spectrum  of  a musical  selection  is  not  only 
non-zero,  but  not  even  flat. 

Turning  our  attention  to  (5.5b),  the  problem  is  not  whether  the 
last  term  averages  to  zero  but  rather  computing  the  average  in  a 
meaningful  way.  The  problem  is  in  deciding  the  actual  value  of  the 
phase.  Values  of  the  arctangent  function  used  in  this  computation 
are  between  0 and  2n;  the  actual  value  may  differ  from  this  by  any 
integer  multiple  of  2n.  The  problem  of  "unwrapping  the  phase"  or 
computing  its  proper  value  is  complex  and  the  object  of  current 
research. 

Ue  will  concentrate  on  the  average  magnitudes  of  (5.5a). 
Experimentation  has  shown  that  the  ear  is  relatively  insensitive  to 
phase  (35).  Consequently,  a system  restored  using  an  estimate  of 
only  the  magnitude  of  the  degrading  system  will  nearly  always  be 
subjectively  the  same  as  a more  accurate  system  computed  using  phase 
inf ormation. 

Returning  to  (5.5a),  we  note  specifically  that 
log!Si(f)|  - ( 1/2)  log  ISk(f ) I*  - ( 1/2)  1 og  I s ( f ) where  Is(f)  is  the 


: 


L I 

J ' 

f / , . 
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periodogram  associated  with  £(t).  Thus  (i/n)  log  IS; ( f ) I = (i/2)Ls(f), 
the  averge  log  spectrum  of  s(t).  Similarly,  (i/n)  ^^log  IP,(f ) I is  the 
ALS  of  p ( t ) and  (5.5a)  can  be  written 

(l/2)LP(f)  * I03IH  ( f ) I + <U2)Ls(f)  . (5.  G) 

To  compute  an  estimate  of  log|H(f)|,  we  make  the  following 
assumption.  For  a modern  recording,  the  first  term  on  the  right 
side  of  the  equal  sign  in  (5.6)  is  a constant  since  for  practical 
purposes  the  recording  equipment  has  a flat  response.  I f we  assume 
tnat  the  spectrum  of  a modern  recording  is  similar  to  that  of  the 
acoustic  recording,  then  an  average  such  as  (5.6)  for  the  modern 
recording  can  be  subtracted  from  (5.6)  leaving  the  first  term  alone. 

If  LM(f)  is  the  average  log  spectrum  of  a modern  prototype 
recording,  then  we  will  define  an  estimate  of  loglH(f)l  as 
log  IH ' ( f ) I = d/2)  LP  ( f ) - (i/2)L„(f) 

a logIH  ( f ) 1 + (1/2)  Ls  ( f ) - ( 1/2)  L„ ( f ) (5.7) 

where  1 og  IH ' ( f ) I is  an  estimate  of  loglH(f)|.  Applying  the  results 
of  (3.12),  we  can  easily  compute  the  expectation  and  variance  of 
this  estimate  as 

EdogIH' (fill  « logIH ( f ) I 

+ (l/2)  logGs  ( f > - (1/2)  logGM ( f ) >=  logIH  ( f ) I (5.8a) 

and 

var  (logIH' (f ) |)  2 var  { (1/2) Ls  ( f ) ) + var  ( ( 1/2) LM  ( f ) ) 

- 2-(itJ/24N)  = nJ/12N  - 0.82247---/N  . (5.8b) 

5.4  Power  Spectrum  Deconvolution 

The  above  approach  suggests  that  an  alternative  estimate  could 
be  obtained  using  the  log  average  spectrum  rather  than  the  averge 
log  spectrum.  Doing  this  (5.7)  becomes 


i 

; 


i 

i 
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log  |H ' ( f ) I = d/2)  pp(  f ) - d/2)p„(f) 

- log  |H  ( f ) I + { 1/2)  p s ( f ) - ( 1/2)  pit  ( f ) . (5.9) 

The  expected  value  and  variance  for  this  approach  are 

ElloglH'(f)l)  a log  |H  ( f ) I + Gs(f)  - GM(f) 

= log  |H  ( f ) I (5.10a) 

and 

var  (log  IH  (f ) |!  a varips(f)l  + var(pM(f)f 

- 2*ty'(N.'/4]  a 1/2N  = 0.5/N  . (5.10b) 

From  (5.9)  and  (5.10)  we  see  that,  for  the  stationary  noiseless 

case,  both  approaches  yield  unbiased  estimators  of  log  IH  ( f ) I . 
However,  the  log  average  estimate  is  more  stable. 

The  first  method,  involving  the  ALS  estimator,  is  refered  to  as 
the  homomorphic  estimator  while  the  second  is  called  the  power 
spectrum  estimator.  Clearly,  except  for  their  variance,  the  methods 
are  equivalent  for  stationary,  noise-free  signals.  In  both  cases, 
the  restoration  filter  will  be  the  inverse  of  |H(f)|. 

5.5  The  Effect  of  Additive  Noise 

In  deriving  (5.8)  and  (5.10),  we  have  neglected  the  effects  of 
noise  in  the  deconvolution  procedure.  Returning  to  (5.1),  and  now 
considering  the  noise,  n(t),  we  again  window  and  compute  the  Fourier 
transform  and  complex  logarithm  giving 

logPi(f)  a log  [S,(  f ) -H  ( f ) + N ( f ) 3 . (5.11) 

Because  of  the  sum  under  the  brackets  in  (5.11),  the  right-hand  side 
of  the  equation  does  not  reduce  to  a sum  of  logarithms.  Proceeding 
anyway,  we  again  define  an  estimate  of  loglH(f)|  as  in  (5.7)  and 
(5.9)  so  that  for  the  homomorphic  approach  we  have 
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log  IH'  ( f ) I = ( 1/2)  LP  - (i/2)LM(f)  (5.12) 

and  for  the  power  spectrum  approach 

log  IH'  ( f ) I = (i/2)pp(f)  - d/2)p„(f)  . (5.13) 

Computing  the  expectation  gives  us 

E I log  IH'  ( f ) I)  s (1/2)  logGP(f)  - (i/2)  logG„(f)  (5.14) 

for  both  approaches.  Now,  if  the  noise  is  uncorrelated  from  the 
signal,  v(t),  then  Gp(f)  = Gs(f)-IH(f)  I2  + G„(f)  and  (5.14)  becomes 
E (log  IH'  (f ) II  a 

a d/2)  log  t (Gs(f)-IH(f ) I*  + GN  ( f ) ) /Gs(f)]  (5.15) 

where  we  have  again  assumed  that  the  spectrum  of  the  modern 
prototype,  G„(f),  equals  the  spectrum  of  the  original  signal,  Gs(f). 

From  (5.15)  we  are  thus  motivated  to  form  the  compensating 
f i 1 ter,  R (f ) , as 

R ( f ) - exp  [(-i/2)  log  [ (G;,  ( f)  -|H(f)  I2  + G(J  ( f ) ) /Gs  ( f ) ) > 

= (Gs  ( f ) / CGS  (f)-IH(f)|2  + G„(f ) 1 ) 1/2  . (5.1G) 

He  see  that  for  the  noise-free  case,  (5.16)  reduces  to  the  inverse 
of  IH(f)|.  In  the  noisy  case,  (5.16)  has  the  interesting  property 
of  naturally  preventing  ill-conditioning,  i.e.,  attempting  to 
restore  a signal  at  frequencies  where  noise  predominates  thereby 
amplifying  the  noise.  However,  the  compensating  filter  of  (5.16) 

becomes  small  when  the  spectrum  of  the  noise  is  large  thus 
preventing  ill-conditioning. 

5.6  The  Effect  of  Nonstctionari ty 

If  the  process  is  now  assumed  to  be  nonstationary,  a much  more 
realistic  assumption,  then  the  above  results  are  modified  in 
accordance  uith  the  results  of  chapter  4.  At  this  point,  assuming 
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that  both  the  modern  anc  acoustic  recordings  are  nonstationary,  we 
model  the  acoustic  recording  as  having  been  passed  through  a 
time-varying  linear  system  with  frequency  responses  (Mf)  and  the 
modern  recording  through  a system  with  frequency  responses  av(f). 

If  the  prototype  recording  is  the  same  selection,  and  the 
recordings  are  noise  free,  then  it  is  not  unreasonable  to  assume  the 
/3t(f)  = a,(f)  for  i = 1,  2,  3,  •••,  N.  In  this  case,  (5.8)  and 
(5.10)  become 

ElloglH'(f)l)  * log  |H(f ) I (5.17) 

for  both  approaches,  but 

var  1 logIH ' ( f ) II  as  n2/12N  (5.18) 

for  the  homomorphic  approach,  and 
var  (log  |H'  ( f ) ) ! as 

- (1/2)  ([(i/N)ZiNillMf)H  / [N((i/n)  I,”  l&mi1)2])  (5.19) 

for  the  power  spectrum  approach. 

Thus,  even  for  nons  tationary  processes,  the  two  approaches  are 
unbiased.  However,  the  power  spectrum  estimate  may  now  be  less 

stable. 

In  general,  however,  the  acoustic  recording  is  not  noise-free 
while  the  modern  recording  is  to  a reasonable  approximation.  Thus, 
absorbing  the  noise,  as  before,  in  the  &(f)  terms,  (5.17)  becomes 
EllogIH'  ( f ) |)  * log  IH  ( f ) I 

+ (1/2)  C (1/n)  Z^logl^m  I2  - (1/N)  2&l«itfH*J  (5.20) 
for  the  homomorphic  estinator,  and 

E ( log  IH'  ( f ) II  * d/2)  log  IH  ( f ) I 

+ (1/2)  Clog  (l/N)  2iNilk(f)  |2  - logd/N)  2,\la,(f ) I2)  (5.21) 

for  the  power  spectrum  estimator.  As  discussed  in  chapter  5,  the 


91 


homomorphic  estimator  tends  to  be  influenced  more  by  the  presence  of 
additive,  stationary  noise.  Thus,  we  would  expect  the  estimator  in 
(5.20)  to  exhibit  more  bias.  As  shown  in  section  5.7,  this 
observation  is  dramatically  observed  in  actual  computations. 

5.7  Experimental  Results 

Two  restorations  are  discussed  in  this  section.  The  first  is  a 
restoration  of  the  190'?  recording  of  Caruso  singing  "Vesti  la 
Giubba"  depicted  in  figure  4.1.  The  other  is  the  simulated 
nonstationary  signal  from  figure  4.7.  This  latter  restoration  has 
the  property  that  the  original  resonating  system,  H(f),  is  availble 
for  direct  comparison  with  the  estimates. 

For  both  these  experiments,  the  data  from  figure  4.2  (Jussi 
Bjoerling  singing  "Vesti  la  Giubba")  was  used  as  the  modern 
prototype.  As  can  be  seen  from  the  activation  spectrum,  figure 
4.2(c),  this  signal  is  relatively  noise  free.  Before  using  the  log 
spectral  estimates  for  this  recording,  however,  they  were  further 
smoothed  in  frequencies  by  convolving  the  estimators  with  a spectral 
window  (as  discussed  in  section  4.G).  This  is  justified  since  the 
additional  smoothing  reduces  the  variance  of  the  system  estimate  and 
no  major  resonance  phenomena  are  expected  in  this  data  which  could 
produce  sharp  peaks.  Figures  5.2(a)  and  (c)  show  both  the  ALS  and 
LAS  estimates.  Figures  5.2(b)  and  (d)  are  these  estimates  smoothed 
in  frequency. 

Figure  5.3  shows  the  esimates  of  the  resonant  frequency 
response,  1 og IH ( f ) 1 , and  the  compensating  filters  from  (5. IB)  as 
computed  by  both  approaches.  Note  that  these  restoration  filters 
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are  truncated  outside  the  effective  bandwidth  of  the  process.  This 
is  to  further  prevent  ill-conditioning  and  to  reduce  the  surface 
noise  of  the  recording  as  much  as  possible. 

There  is  clearly  a difference  in  the  two  compensating  filters. 
As  predicted  by  (S.20)  ard  (5.21),  this  is  expected.  Because  of  the 
sharp  resonant  peak  in  the  log  spectral  estimates  of  the  Caruso 
recording,  the  biasing  effect  of  the  surface  noise  affects  the 
homomorphic  estimator  principly  in  the  low  and  high  frequency 
regions  with  the  effect  being  most  pronounced  for  the  bass 
frequencies.  The  failure  of  the  homomorphic  filter  to  properly 
compensate  in  the  bass  region  is  apparent  in  figure  5.3  (b). 
However,  as  there  is  nothing  to  compare  this  filter  with,  it  is  not 
possible  to  measure  this  bias  and  verify  that,  in  fact,  the 
homomorphic  estimate  is  biased  more.  It  is  also  clear  that  the 
homomorphic  estimate  is  more  stable.  Again,  this  is  predicted  by 
(5.18)  and  (5.19). 

To  compensate  for  the  bias  in  the  bass  region  of  the 
homomorphic  filter,  an  empirical  bass  boost  was  added  (figure  5.4). 
As  is  evident  from  figure  5.4(b),  this  has  the  effect  of  producing  a 
filter  much  closer  to  the  power  spectral  filter,  figure  5.3(d). 

Auditioning  of  restorations  produced  by  these  two  filters  also 
reveals  a difference.  Unquestionably,  both  restorations  show  a 
definite  improvement  in  the  resonant  quality;  the  reverberations  so 
obvious  in  the  original  recording  are  missing.  In  general,  however, 
the  power  spectral  restoration  is  more  pleasing  to  the  ear. 
Interestingly,  though,  not  all  restorations  attempted  by  these  two 
methods  exhibit  the  same  preferential  ordering.  It  is  apparent  that 
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the  presence  of  noise  effects  the  listening  quality  of  the 
restorations  in  a more  complex  way  than  just  biasing  the 
compensating  filter. 

Figures  5.5  and  5.5  are  the  log  spectral  estimates  for  the 
homomorphic  and  power  spectrum  restorations,  respectively.  It  can 
be  seen  that  the  effect  of  filtering  the  Caruso  signal  is  to  impart 
to  it  the  general  shape  of  the  spectrum  of  the  modern  prototype 
(compare  to  figure  A. 2).  The  effect  of  the  nonstationary  biases  is 
also  evident  in  the  fac:  that  the  ALS  estimate  of  the  homomorphic 
restoration  is  very  smooth  compared  to  the  LAS  estimate.  This 
results  from  the  details  of  the  ALS  estimate  used  to  produce  the 
compensating  filter  cancelling  when  computing  this  ALS  estimate; 
similar  results  occur  for  the  power  spectrum  restoration  and  the  LAS 
estimates. 

Figure  5.7  shows  the  results  of  a similar  restoration  for  the 
simulated  process  of  figure  4.7.  For  this  data,  however,  the 
resonant  system  is  known  (see  figure  4.7(b))  and  can  be  compared  to 
the  restoration  filters,  figure  5.7(b)  and  figure  5.7(d).  Figure 
5.8  is  the  sum  of  the  compensating  filters  and  the  resonant  system 
and  represents  the  bias  in  the  filters.  The  bias  in  the  homomorphic 
filter,  figure  5.8(a)  is  obvious.  The  fact  that  the  two  simulated 
filters  in  figures  5.7(b)  and  (d)  are  similar  to  the  actual  filters 
for  the  Caruso  recording  is  further  support  of  the  effectiveness  of 
the  model  we  adopted  of  nonstationarity,  and  strongly  supports  the 
conclusion  that  the  homomorphic  approach  is  influenced  more  by 
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FIGURE  5.2 

Log  spectral  estimates  (N  = 544)  for  the  modern  prototype  from 
figure  4.2  used  in  estimating  the  resonant  system,  H(f):  (a)  average 
log  spectrum,  (b)  average  log  spectrum  smoothed  in  frequencies,  (c)  log 
average  spectrum,  and  (d)  log  average  spectrum  smoothed  in  frequencies. 


FIGURE  5.3 

Estimate  of  log|H(f)|  for  the  Caruso  recording  of  figure  4.1:  (a) 
homomorphic  estimate  of  log |H( f )|  , (b)  homomorphic  compensating  filter 
(c)  power  spectrum  estimate  of  log|H(f)|,  and  (d)  power  spectrum  com- 
pensatina  filter. 
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FIGURE  5.4 

Emperical  bass  boost  to  compensate  for  noise  biasing  the  homo- 
morphic compensating  filter,  figure  5.3(b):  (a)  bass  boost  and  (b) 
homomorphic  compensating  filter  with  the  emperical  bass  boost  added. 
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FIGURE  5.5 

Log  spectral  estimates  (N  = 470)  for  the  homomorphically  restored 
Caruso  recording,  figure  4.1:  (a)  log  average  spectrum,  (b)  average  log 
spectrum,  (c)  activation  spectrum,  and  (d)  spectral  dynamic  range. 
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FIGURE  5.7 

Estimates  of  log[H(f)j  for  the  simulated  nonstationary  process  of 
figure  4.7.  Compare  with  the  actual  resonant  system,  figure  4.8(b): 
(a)  homomorphic  estimate,  (b)  homorphic  compensating  filter,  (c)  power 
spectrum  estimate,  and  (d)  power  spectrum  compensating  filter. 
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FIGURE  5.8 

Bias  in  the  compensating  filters  of  figure  5.7  computed  as  the  sum 
of  the  resonant  filter,  figure  4. 8(b)  and  the  filters  in  figure  5.7(b) 
and  figure  5.7(d):  (a)  homomorphic  bias  and  (b)  power  spectrum  bias. 
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CHAPTER  e 


CONCLUSION  AND  SUMMARY 


(G. lb) 
(G.lc) 
(G. Id) 


G.l  Statistical  Summary 

The  principle  equations  representing  the  statistical  properties 
of  the  log  average  spectrum  (Lx(f))  and  average  log  spectrum  (px(f)) 
are 

E { ( f ) 1 at  logGx(f)  + ^(N)  - log(N)  * logG*(f),  (G.la) 
var(Mf))  a <P'  IN)  a 1 / N, 

EfLx(f)I  at  logGx(f)  - y, 
var!Lx(f))  a ti1  / (GN)  a:  l.G443>*»/N, 

and 

E(Ax(f))  = log(N)  - •// (N ) - y a -y 
for  the  stationary  process  x(t)  and 

E (Lx ( f ) ) at  iogG^If)  - y + (i/n)  2i"log|&(f)  I2, 
var  ILX ( f ) ) at  nVGN  at  1.6449—/N, 

E IPx(f ) ) a:  logG„(f)  + log  E(i/n)  ZAI0.(f)  1*3 , 
var(Px(f)i  at 

t (1/N)  2Aia.(f)l43  / (N*  [ (l/N)  X&l0i<f)|*3*)t 


(G.le) 


(6.2a) 
(G. 2b) 
(G. 2c) 


(G. 2d) 


and 


E (Ax( f ) ) a -y 

+ ( (i/N)  Iogl0t(f)|j  - log  [ (l/N)  JAUMf)  1*3 ) (G.2e) 

for  a nonstationary  process  x(t)  derived  from  the  stationary  process 
y(t).  In  all  equations,  Gx(f)  and  Gw ( f ) are  the  spectral  density 
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functions  of  x ( t ) and  y(t),  respectively  and  (3t(f)  is  the  frequency 
response  of  the  ith  time-varying  linear  system  (see  section  4.2). 

6.2  Practical  Conclusions 

The  results  of  this  research  have  practical  significance.  The 
log  average  spectrum  commonly  arises  in  conventional  spectral 
analysis;  the  average  log  spectrum  arises  naturally  in  certain 
applications  of  homomorphic  signal  processing  and  is  an  interesting 
alternative  for  estimating  the  log  spectrum  of  a signal.  Clearly, 
each  estimator  has  advantages  (and  disadvantages)  that  should  be 
considered  for  any  particular  application. 

The  principal  advantages  of  the  log  average  spectrum  are  (1)  a 
faster  computation  time  (18%  less  than  the  average  log  spectrum  on  a 
PDP-10  computer  system),  (2)  it  is  a smoother  estimate  of  stationary 
processes,  an^l  (3)  for  r.onstationary  processes  it  is  affected  less 
significantly  by  additive,  stationary  noise.  Its  disadvantages  are 
(1)  for  nonstationary  processes,  it  tends  to  be  less  stable  and  (2) 
it  has  a lower  coherent  signal  to  noise  ratio. 

Similarly,  advantages  of  the  average  log  spectrum  are  (1) 
generally,  it  will  exhibit  more  stability  for  nonstationary 
processes  and  (2)  it  has  a higher  coherent  signal  to  noise  ratio. 
Disadvantages  are  (1)  it  has  less  stability  for  stationary 
processes,  (2)  it  can  be  significantly  affected  by  the  presence  of 
additive  noise  for  nonsi ationary  processes,  and  (3)  it  requires  a 
longer  computation  time. 

It  is  the  suggestion  of  the  author  that  in  general  both 
estimators  be  computed  in  practical  spectral 


analysis. 


The 


103 


additional  computation  time  is  not  significant  considering  the 
relative  advantages.  By  doing  this,  not  only  can  the  two  estimators 
be  compared,  but  the  activation  spectrum  can  be  computed.  As  shown 
in  chapter  4,  this  can  reveal  significant  and  interesting 
information  about  the  nature  of  a process. 

In  the  situation  where  these  estimators  are  used  to  estimate  a 
system  response,  as  in  chapter  5,  the  LAS  will  generally  give  better 
results  since  it  is  least  affected  by  noise.  However,  in  any 
particular  experiment,  either  approach  is  potentially  the  more 
desirable  in  terms  of  achieving  the  desired  goals. 

In  actually  computing  the  LAS  and  ALS,  two  considerations  are 
important.  First,  the  choice  of  a spectral  window  affecting  the 
overall  bias  of  the  estimates.  It  is  important  to  choose  a window 
that  will  yield  the  desired  resolution  (narrow  peak)  yet  bias  the 
estimates  as  little  as  possible  (small  side  lobes).  Unfortunately, 
these  two  criteria  often  conflict  (271. 

Second,  the  choice  of  the  data  segment  lengths.  Frequently, 
this  choice  is  constrained  by  minimum  resolution  requirements  and 
the  total  amount  of  data  available.  Generally,  it  is  best  to 
compute  the  minimum  length  required  to  give  the  desired  resolution 
and  use  this  to  determine  the  number  of  segments  to  be  used  in  the 
smoothing  process.  Keep  in  mind  that  the  variance  of  the  estimators 
can  often  be  further  reduced  by  overlapping  the  data  segments  as 
proposed  by  Welch  [22). 

G.3  Further  Research 

There  are  several  areas  of  additional  research  of  both  a 
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theoretical  and  practical  nature.  The  results  derived  in  this 
research  are  for  spectral  estimators  smoothed  by  the  Bartlett 
averaging  procedure.  While  there  is  reason  to  believe  similar 
results  apply  to  estimates  smoothed  by  other  techniques,  it  uould  be 
useful  to  extend  this  analysis.  For  example,  one  might  consider  the 
effect  of  smoothing  the  log  periodogram  by  convolution  with  a 
spectral  window.  Similarly,  the  results  could  be  extended  to 
include  the  case  where  the  data  segments  overlap. 

Finally,  it  would  be  desirable  to  extend  this  analysis  into  two 
dimensional  signal  processing  (images).  In  fact,  these  issues  have 
been  encountered  in  some  image  processing  research.  For  example,  in 
his  doctoral  dissertation  [24] , Cannon  computes  both  the  LAS  and  ALS 
estimators  of  an  image.  As  shown  in  figure  G.l  (reprinted  from  [24] 
with  permission  from  T.  H.  Cannon)  there  is  clearly  a difference 
between  the  two  log  spectral  estimators  (note  that  the  activation 
spectrum  was  not  explicitly  computed).  It  is  reasonable  to  believe 
that  the  results  developed  for  one  dimensional  processes  generally 
apply  in  two  dimensions  since  the  mathematics  involved  can  be 
readily  extended  into  two  dimensions.  Computation  of  the  activation 
spectrum  for  an  image  will  undoubtedly  enhance  the  understanding  of 
the  image  as  well  as  orovide  useful  insight  into  selection  of 
prototypes  for  image  deblurring  as  proposed  by  Cannon  [24]  , 
Cole  (253  and  others. 
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FIGURE  6.1 

Log  spectral  estimates  for  an  image  (reprinted  with  permission 
from  [24]):  (a)  the  image,  (b)  log  average  spectrum,  and  (c)  average 
log  spectrum. 
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APPENDIX  A 
SPECIAL  FUNCTIONS 

A.l  Euler’s  Constant 

The  Euler-Mascheroni  constant,  Y,  or  simply  Euler’s  constant, 
is  given  by  137, p. 9] 

y - tiff  A„  (A. la) 

where 

Am  - Zui  (i“‘)  - log  (N)  . (A.  lb) 

Numerically,  Y - 0.57721  56649  01533—  . 

A. 2 The  Gamma  Function 

The  gamma  function,  T(t),  is  defined  as  [37, p. 8) 

T(t)  = J'"xt‘‘*6!xp  (-x) dx  . (A.  2) 

By  partial  integration  of  (A. 2),  we  have 

r(t+l)  - t-r(t)  (A. 3) 

and  by  substitution 

r (i ) - r(2)  - l . (a. 4) 

A. 3 The  Digamma  and  Tricamma  Functions 

The  so-called  digarnma  function  (or  Euler’s  psi  function),  ^(t), 
is  the  logarithmic  derivative  of  the  gamma  function  [37, p. 121, 

yMt)  - (d/dt)  logT(t)  - r'(t)/r(t)  . (A. 5) 

There  are  several  representations  of  ^(t)  [37, pp. 12, 13, 16)  . Among 
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\Mt+l)  - -y  + 2i"  d’1  - (t+i)'*)  . (A.B) 

For  N - 1,  2,  3,  — , we  can  write 

yMN+t)  = 1/t  + l/(t+l)  + — + 1 / (t+N-1)  + \Ht)  . (A. 7) 

From  (A. G) , we  have 

*(1)  - -v  - -0.57721—  (A. 8) 

giving 

yMN+1)  - 1 + 1/2  + 1/3  + — + 1/N  - v . (A. 9) 

Approximations  from  asymptotic  expressions  for  \Mt)  are  useful. 
Using  a Euler-tlaclaurin  expansion  [33, p. 483) 

iMt)  = log ( t ) - 1 / (2 1)  - 1 / ( 1 2 1 2 ) + 

1/  <120t4)  - •••  (A. 10a) 

so  that 

yMt)  ^ log(t)  as  t -*  oo  . (A. 10b) 

A more  accurate  approximation  is  given  by  Cox  and  Lewis  [20,p.26)  as 

*Mt)  - log  ( t ) - 1 / (2 1 - 1/3  + 1 / I16t] ) • (A. 11) 

However,  (A. 10b)  is  suitable  for  our  needs. 

The  first  derivative  of  the  digamma  function,  ^'(t),  is 
frequently  encountered  and  called  the  trigamma  function.  By 
di  f ferentiating  (A.  6),  we  have  [37,p.26) 

V^(t+1)  =*  SWt+i)'2  (A.  12a) 

and 

^'(1)  - X"i<i>*2  - it2 / 6 - 1.G449  34067—  . (A. 12b) 

By  differentiating  (A. 7),  for  N = 1,  2,  3,  — 

vfrMN+l)  - rt2/6  - 1 - 1/4  - 1/9  - — - 1/N2  . (A. 13) 

As  an  approximation  for  ^'(t),  from  (A. 10a)  we  derive 
^'(t)  - (1  / t ) • Cl  + 1 / ( 2 1 ) + 


1/  (6t2)  - 1/  ( 30 1 * ) 


— ] . 


(A. 14) 
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from  which  we  conclude 

^'(t)  * 1/t  as  t - ® , , (A. 15) 

Bartlett  and  Kendall  C333  give  a sightly  more  accurate  approximation 

Vfr'(t)  * 1/  (t  - 1/2)  (A.1G) 

and  Cox  and  Lewis  [20,p.26I  give 

Vfr'(t)  * 1 / ( t - 1/2  + 1 / llOtl  > . (A. 17) 

Table  A.l  is  a list  of  values  of  \MN) , log(N),  ^'(N),  and  1/N 
for  selected  values  of  N.  Comparison  of  appropriate  values 

indicates  the  asymptotic  accuracy  of  (A. 10b)  and  (A. 15).  Figure  A.l 
is  a graphical  representation  of  the  data  in  Table  A.l.  See 

[AO, pp. 945-61  for  additional  integral  and  series  representations  of 

rit),  iMt)  and  ( t ) . 


TABLE  A.l 


SELECTED  VALUES  OF  SPECIAL  FUNCTIONS 


N 

iMN) 

Log  (N) 

\f/'  (N) 

1/N 

1 

-0.577216 

0.000000 

1.644934 

1.000000 

2 

0.422784 

0.G93147 

0.644934 

0.500000 

3 

0.922784 

1.098612 

0.394934 

0.333333 

4 

1 . 256118 

1.38G294 

0.283823 

0.250000 

5 

1 .5061 18 

1.609438 

0.221323 

0.200000 

G 

1.706118 

1.791759 

0.181323 

0.1666G7 

7 

1.872784 

1.945910 

0.153545 

0.142857 

8 

2.015G41 

2.079442 

0.133137 

0.125000 

9 

2.140G41 

2.197225 

0.117512 

0.111111 

10 

2.251753 

2.302535 

0.105166 

0.100000 

20 

2.970524 

2.995732 

0.054041 

0.050000 

30 

3.384438 

3.401197 

0.033895 

0.033333 

40 

3.G7G327 

3.688879 

0.025315 

0.025000 

50 

3.901990 

3.912023 

0.020201 

0.020000 

100 

4.G001G2 

4.605170 

0.010152 

0.010000 

APPENDIX  B 


DERIVATIONS 


B.l  Log  Chi-square  Statistics 

The  log  chi-square  cistribution  was  first  described  by  Bartlett 
and  Kendall  [391  and  has  subsequently  been  used  by  others  (e.g.,  see 
[133,  1203,  1253,  and  [303).  It  has  particular  application  in  the 

statistical  analysis  of  the  logarithm  of  spectral  estimators. 

The  probability  density  function  of  a log  chi-square  random 
variable  may  be  easily  derived  from  the  chi-square  dis tribut ion. 
Let  Y - g (X)  =■  log(X)  where  X = r*x2n  and  X2n  is  a random  variable 
having  a chi-square  distribution  with  n degrees  of  freedom.  Then 
Y - log(X)  =»  log(r'X7n)  has  a log  chi-square  distribution. 

The  pdf  of  X is  easily  obtained  from  (2.11)  and  is 
f„(x)  = [r-2n,2r  (n/2)  3 '*• 

(x  / r)  n/2"‘-exp  (-x  / [2r3  ) . (B.l) 

To  determine  the  pdf  cf  Y,  we  note  that  X = g~‘(Y)  = exp(Y)  = 

(d  / dy) g‘l  (y) . Equation  (2.G)  may  be  used  so  that 
f„(  y)  = Mg"‘(y)H(d/dy)gM(y)l 

- [r»2T  (n  / 2)  3 **•  lexp  (y)  / (2r)  3 


exp (-exp  (y)  / [2r] ) -exp (y) 

= r (n  / 2)  "*•  [exp  (y  - log  (2r) ) 3 

exp(y  - log(2r)  - exp  [y  - log(2r)3)  (B.2) 

where  the  last  step  is  obtained  by  noting  that  (2r)"‘  = 


Ill 

exp C-log (2r) ] . For  n - 2,  the  first  two  factors  are  unity  and  (B.2) 
reduces  to 

f s (y ) = exp(y  - log(2r)  - exp [y  - log(2r)])  . (B.3) 

If  we  note  that  E IX)  = //x  = Elr^x^l  = r-n  (see  Table  2.1),  then  for 
n » 2,  Mx  = 2r  and  (B.3)  becomes 

f y ( y ) “ exp  (y  - log//*  - exp  [y  - log//xl)  (B.4) 

as  given  in  (2.13)  . (13.4)  has  been  derived  by  others  [411  , 125) 

directly  from  an  exponential  distribution. 

The  mean  and  variance  of  V = log(r-x2r)  a<~e  derived  by  Bartlett 
and  Kendall  (391  from  I1a(t),  the  characterist  ic  function  of  V. 
fly(t)  is  given  by 

f1u(t)  = Elexp(jtY))  = Elexp(  jtlog(X)  ))  = E (XJt) 

■=  J':xjt(r-2r(n/2)]'l-(x/  [2rl  ) n/J-‘-exp  (-x  / (2rl ) d* 

= r(n/2)'M2r)  jtJ':zJ,*n/l-‘-exp(-z)dz  (B.5) 

where  z - x/(2r)  in  the  last  step.  Noting  from  (A. 2)  that 

r (n  / 2 + jt)  = j'»zjt*n/j-i<exp  (z)  dz,  (B.5)  then  becomes 

riy(t)  - (2r) [T  (n  / 2 + j t)  / T (n  / 2)  3 <B . G) 

and  thus  the  cummulant  function  of  Y is 
Ky(t)  = log  (nit)) 

- j t*  log  (2r)  + log  T (n  / 2 + jt)  - logT(N/2)  . (B.7) 

Hence  the  mean  and  variarce  are 

E (Y)  - j*‘(d/dt)Ky(t)  I (t»0)  ■=  log  (2r)  + \Mn/2)  (B.8a) 

and 

variY)  - j"2  (d2/  dtl)  Ku  ( t ) I ( t=0)  = *'(n/2>  (B.8b) 

where  ^(t)  and  ^'(t)  are  the  digamma  and  trigamma  functions 
respectively  (see  section  A. 3).  For  n = 2,  and  again  noting  that 


E (X)  - Mx  - 2 r,  (B.8)  can  be  written 
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E IV)  = log(2r)  + y''(l)  = log//*  - V 

and 

var(Y)  - ^'(1)  - n'/B  - 1.6449...  . 


(B. 9a) 
(B. 9b) 


B.2  The  Periodogram 

Let  x(t)  be  a zero-mean,  stationary  random  process  and  xT(t)  a 
sample  function  of  x(t)  cn  the  interval  [0,T).  Then  the  periodogram 
of  x(t)  is  defined  as 

I*(f)  = flcx^tr))  = J^CxX(t) -exp  {-2n  j f r)  dr  (B.10) 

where  Cxx(t)  is  the  sample  autocovariance  function  given  by  (2.33a). 
An  alternate  definition  of  Ix(f)  in  terms  of  the  Fourier  transform 
of  xT(t)  may  be  derived  as  follows. 

Write  (B.10)  as 

Ix(f)  - x?  c Cl/T)  J'rT,Xr(t).xr(t+lr|)dt].exp(-2njfT)dT 

+ XJ  C (l/T)  J':-'r,xT(t).x(t+lr|)dt).exp(-2njfr)dT  (B.ll) 

where  (2.33a)  has  been  substituted  for  cXx(t)  with  fi*  = 0.  By  using 
the  transformation  of  variables 

s=t+r  (B. 12) 

and  appropriately  rewriting  the  limits  of  integration,  (B.ll) 
becomes 


Ix(f)  = (l/T)  • J’ J"oXT  ( t ) .xT  ( s) -exp  (-2n  j f (s- tl  ) dt  ds 
. = (i/t)  • X;  xT  (s) -exp  (-2n  j f s)  ds- 
Xo><"  ( t)  - exp  (+2rt  j f t ) dt 

= (i/T)  -XT  ( f ) -XT  ( f ) * = ( i/T ) - |XT  ( f ) I2  (B.  13) 

where  * denotes  complex  conjugation  and  XT(f)  is  the  finite  Fourier 
transform  of  xT(t).  (B.13)  is  the  same  expression  as  given  in 
(2.41).  A more  detailed  derivation  is  given  in  [4, pp. 214-15)  and 
C3,pp. 82-84] . 


113 


The  definition  of  Ix(f)  in  (B.13)  can  be  used  to  derive  the 
distribution  of  the  periodogram.  For  a zero-mean,  Gaussian,  purely 
random  (white)  process,  this  can  be  done  precisely  and  rather 
simply.  Results  for  other  processes  are  mathematically  complex,  but 
have  been  derived  (e.g.,  see  (42) ) . For  correlated  processes, 
however,  these  simpler  results  hold  asymptotically  and  even  for 
non-Gaussian  processes  are  generally  quite  accurate. 

Let  x ( t ) be  defined  as  before  with  the  further  restriction  that 
it  be  Gaussian  and  white.  Define 

XT(f)  = J'oXr  ( t ) *exp  (-2n  j f t ) dt 

= Xo *T  ( t ) -cos  (-2nf  t ) dt  + ]•  Xo*T  ( t ) ‘sin  (-2nf  t ) d t 
- X„(f)  + j.X,(f)  (B.  14) 

where  X„(f)  and  X,(f)  denote  the  real  and  imaginary  parts  of  XT(f) 
respectively.  Since  xT(t)  is  normal  for  each  t,  0 < t < T,  and 
linear  combinations  of  weighted  normal  random  variables  are 
themselves  normal,  it  follows  that  X„(f)  and  X,(f)  have  Gaussian 
distributions.  For  a non-Gaussian  process,  this  will  be 

approximately  true  by  the  Central  Limit  Theorem.  Furthermore, 

E (X„  ( f ) ) - E IX,  ( f ) 1 = 0 . 

Koopmans  [8, pp. 2G1-E3) , Jenkins  and  Watts  [4,p.239]  and  others 
show  that  at  the  harmonic  frequencies  f = k/T,  |k|  = 0,  1,  2,  ••• 

(and  at  all  frequencies  for  large  T)  , the  vector  pairs 

(X,  ( f ,)  , X,  ( f ,) ) and  (X,  ( f?)  , Xx  ( f ,) ) are  independent  for  f,  * f?. 
This  follows  from  the  fact  that 

cov  (X,(f,)  ,X„(f?) ) - covlX,(f,)  ,X,(fj))  -0  (B.1G) 

and  that  X„(f)  and  X,(f)  are  Gaussian.  It  can  be  similarly  shown 


that  X„(f)  and  X, ( f ) are  independent  of  each  other. 
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Ule  can  now  determine  the  distribution  of  the  periodograni. 
First,  the  variance  of  X„(f)  may  be  easily  derived  for  the  harmonic 
frequencies  (and  again  at  all  frequencies  for  large  T) 

var(X„(f))  = var  f J’ xT  ( t) -cos  (-2nf  t)  dt) 

■ Jo  var  (:<T  ( t) ) «cos2  (2nf  t)  dt 
= <7-V  Jo  cos2  (2nf  t ) d t 

= T ( cr2x  / 2 ) , f = k/T,  Ik  I = 0,  1,  2,  •••  . (B.  17) 

Similarly, 

var  (X,  ( f ) ) = T(<72*/2),  f = k/T,  Ikl  = 1,  2,  •••  . (B.18) 

Now,  from  (B.13)  and  (B.14) 

Ix(  f ) - (i/T)  • [X,2(f)  + X,2(f)]  . (B.  19) 

and  we  see  that  the  periodograni  is  the  sum  of  two  independent, 
identically  distributed,  squared,  normal  random  variables.  Thus  it 
has  a distribution  proportional  to  a chi-square  distribution  with  2 
degrees  of  freedom,  r-x2,.  To  determine  r,  we  note 

E (Ix(f) ) = (1/t)  • (E  IX,2 ( f ) ) + E IX,2 ( f ) ) ) 

= ( i/t)  • ( var  (X,  ( f ) ) + var(X,(f)I) 

= (1/TMT-0-V2  + T-<72,/2)  = <t2h  (B.  20) 

and  using  (2.18b)  and  (2.48) 

r - E (I„(f)  I / n - o'2*  / 2 . (B.21) 

It  f ol lows  then  that 

2-I*(f)  = x2?  . (B.  22) 

It  can  be  shown  (8, pp. 2S5-74]  that  similar  results  apply  to 
smoothed  spectral  estimators  with  the  degrees  of  freedom  a function 
of  the  spectral  window.  Specifically,  the  Bartlett  estimator, 
P,(f)  - (l/N)  is  an  average  of  N x??  random  variables.  Since 


x(t)  is  Gaussian  and  white,  each  xv(t)  is  independent  so  that  each 
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periodogram  is  also  independent  and  thus 

2N-Px(f) = x??N  . (B.  23) 

These  results  are  exact  for  a Gaussian,  white  process.  More 
general  results  are  discussed  by  Jenkins  and  Watts  14], 
Koopmans  [81,  Hannan  [42!!  , and  others,  and  are  much  more  complex. 
However,  the  above  are  asymptotically  valid  for  non-Gaussian, 
non-white  processes  and  (B.22)  and  (B.23)  can  be  rewritten 

2*IX ( f ) /Gx(f)  * x22  (B. 24a) 

and 

2N.px(f)  / Gx ( f ) * X2?H  (B. 24b) 

where  Gx(f)  is  the  spectrum  of  x(t). 

The  equations  in  (B.24)  were  derived  for  the  periodogram  given 
by  (2.40).  However,  they  also  apply  to  Welch’s  modified  periodogram 
(2.58b)  since  weighted  sums  of  normal  random  variables  are  still 
normal . 

Figure  B.l  demonstrates  the  above  discussion  for  a simulated 
Gaussian,  white  process.  Figure  B.l (a)  is  a histogram  for  one  409G 
point  segment  of  the  process.  Superimposed  is  a normal  pdf  computed 
from  the  sample  mean  and  variance  of  the  process.  Although  a 
chi-square  goodness  of  fit  test  (95%  level)  failed  by  a small  margin 
(171.83  compared  to  t lie  test  statistic  120.99),  the  general 
normality  of  the  data  is  apparent. 

Similarly,  Figure  B.l(b)  is  a histogram  of  the  real  part  of  the 
Fourier  transform;  Figure  B.l(c)  of  the  periodogram;  and  Figure 
B.l(d)  of  the  log-periodogram.  Superimposed  are  normal,  X%,  and 
logX?7  distributions,  respectively.  Again,  the  agreement  between 
theory  and  experiment  is  apparent.  In  these  latter  3 figures. 
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chi-square  goodness  of  fit  tests  (95%  level)  showed  agreement  (93.84 
vs.  120.99  for  (b) . 88.67  vs.  122.11  for  (c)  and  107.02  vs. 
122.11  for  d). 


B.3  The  Hyperbolic  Distribution 

I f X is  a random  variable  with  a uniform  distribution  on  the 
interval  Ca,bl,  then  Y » g(x)  ■ exp(X)  has  a hyperbolic 
distribution.  The  form  cf  this  pdf  is  easily  derived  from  (2.6). 

Ue  first  note  that  the  pdf  of  X is  given  by 
f„(x)  «*l/(b-a),  asxsb 

« 0,  otherwise  (B.25) 

Furthermore,  g(X)  has  an  inverse  given  by  X » g‘‘(Y)  - log(Y)  so 
that  | (d/dY)g-‘(Y)  ( - 11 /''I,  Y * 0.  Since  Y =.  exp(X),  however,  Y 2>  0 
and  I (d/dY)g*‘ (Y)  I - 1/Y.  Applying  (2.6),  yields 
f*(y)  - fx(g'‘(y)Hdg'My)  / dy I 

- 1 / ( (b  - a)  y]  , A < y < B 
- 0,  otherwise 
where  A and  B are  given  ty 
A = exp (a) 

and 

B - exp(b)  . 

Before  computing  various  moments  of  Y,  it  is  useful  to  define 
two  quantities  representing  the  dynamic  range  over  which  Y and  X are 
allowed  to  vary.  Specifically 

0 - B/A  (B. 28a) 

and 

d - log (D)  = log (B / A)  » log(B)  - log(A)  «=  b - a (B.28b) 


(B. 26) 


(B. 27a) 


(B. 27b) 


where  D is  the  dynamic  range  of  Y and  d is  its  logarithm. 
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FIGURE  B . 1 

Experimental  histograms  of:  (a)  Gaussian,  white  data  segment;  (b) 
Fourier  transform  (real  part)  of  (a);  (c)  periodogram  computed  from 
(a);  and  (d)  log  periodogram  computed  from  (c).  The  smooth  curves  are 
theoretically  predicted  pdfs  corresponding  to  the  experimental  data. 
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The  expected  value  of  Y is  derived  from  (B.2S)  and  (2.16).  Ue 

have 

E (Y)  - J'Sy(yd)-‘dy  = (B  - A)  / d 

- (B / d)  • (1  - 1/D)  a (B/d)  (B.29) 

where  the  approximation  is  for  large  □.  Similarly,  other  moments 
may  be  computed: 

E(Yl)  = ISyMyd)‘‘dy  = (B*  - A*) /2d 

- (BV2d) • (1  - I/O2)  * BV2d,  (B .30a) 

E (log  (Y) ) - EIX)  - (a  + b)  /2 

- ( 1/2)  log  ( A-B ) =•  { 1/2)  log  (B*  / □) 

■ log(B)  - d/2  = b - d/2,  (B.30b) 

and 

logE(Y)  = log(B/d)  + logtl  - 1/D) 

= log(B/d)  = b - log(d)  (B.3Qc) 


where  again  the  approximations  are  for  large  D. 
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APPENDIX  C 

DETAILS  OF  DIGITAL  INPLENENTATION 


C.l  Digital  Computation  of  the  Periodogram 

The  modified  periodogram  of  (2.58)  may  be  computed  digitally 
for  a discrete  process,  x (n) , by  replacing  the  Fourier  integral  with 
a discrete  Fourier  transform  (DFT) . Equation  (Z.58)  then  becomes 


p„(k)  - (l/N)2AJi(k), 

Jt(k)  = (1  /HU)*l2”'Jx(n)-w(n)*exp(-2njkn/n)  I2, 


and 


(C. la) 
(C. lb) 

(C.lc) 


U - (l/mX”;X<n) 
as  given  by  Uelch  (221. 

The  spectral  window  used  in  this  research  for  (C.lb)  is  one 
form  of  a Tukey  window  called  a Hanning  window  [4, p. 2441.  It  has 

the  continuous  representation 

w(t)  - (1/2)  + (1/2) -cos (nt  /H) , Itl  < M It. 2a) 

in  the  time  domain  and 

U(f)  = M-  (sin  (2nHf ) / 2nMf ) • (1  - [2nf)2)-‘,  If  I < « (C.2b) 

in  the  frequency  domain.  The  discrete  representation  of  (C.2a)  is 

w(n)  - (1/2)  + (1  /2)-cos(2nn/N)  , Ini  < N/2  . (C.3) 

For  this  window,  U - 0.375.  It  is  used  both  in  the  implementation 
of  (C.l)  and  in  the  frequency  smoothing  of  the  log  spectral 
estimators  described  in  section  5.7. 


— a 


The  discrete  Fourier  transform  of  (C.lb)  is  most  efficiently 
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computed  by  using  the  Fast  Fourier  Transform  (FFT)  commonly  used  in 
digital  signal  processing  (e.g.,  see  [21,  [31,  [5],  and  [Gl ) . 

Essentially,  it  makes  uses  of  certain  symmetry  properties  of  the  OFT 
to  reduce  the  computation  time  from  being  proportional  to  N*  to 
N«log?(N).  For  a OFT  of  length  8192  (the  4096  point  segments  in 
this  research  are  augmented  uith  409G  zeros  to  enhance  resolution) 
this  is  a reduction  in  time  by  a factor  of  about  630:1. 

Further  reduction  in  computation  time  is  achieved  by  using  a 
biplexed  FFT.  Using  symmetry  properties  of  the  Fourier  transform  of 
real  data,  this  algorithm  simultaneously  transforms  two  adjacent 
data  segments.  One  segment  is  inserted  as  the  real  part  of  a 
signal,  the  other  as  the  imaginary  part.  After  applying  the  FFT, 
the  individual  transformations  are  then  separated  in  the  frequency 
domain  by  using  odd-even  symmetry.  Time  is  thus  reduced  by  nearly  a 
factor  of  two. 

An  important  consequence  of  digital  implementation  is  the 
effect  of  quantization  and  sampling.  These  topics  are  discussed 

nfbre’~'  fCiTly"  "in’ ~Trttter~  r e~er-emees-  Ae.-g-.  r —see CS14-  . — 2T±ie..-e4  f.ec.t 

campling,  of  course,  is  to  introduce  aliasing  in  the  frequency 
domain.  For  this  reason,  the  data  is  filtered  at  approximately  half 
the  sampling  frequency. 

The  effect  of  quantization  is  to  add  noise.  Uith  a 14-bit  A/D, 
however,  this  is  about  ICO  dB  below  the  signal  level.  Of  course,  in 
the  actual  sampling  process,  it  is  necessary  to  keep  signal  levels 
low  enough  to  prevent  overflowing  the  14-bit  storage  restraint  and 
yet  high  enough  to  maximize  the  number  of  quantization  steps.  In 
some  instances  (such  as  when  filtering  in  the  simulation  of  section 


I 
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5.7),  it  i9  necessary  to  digitally  scale  the  signals  to  prevent 
overflowing  the  finite  length  data  registers. 


i 


C.2  Generation  and  Repre sentation  of  Experimental  Data 

Two  forms  of  data  are  used  in  this  research:  computer  generated 
random  processes  and  real,  physical  data  that  has  been  sampled, 
quantized  and  stored  digitally.  The  computer  generated  data  is  used 
to  test  the  theoretical  nodels  while  the  real  data  represents  actual 
implementation  of  the  theory. 

In  all  cases,  the  data  is  stored  on  high  speed  magnetic  discs 
with  two  14-bit  samples  packed  in  each  3G-bit  word.  The  data  is 
stored  in  integer  formal  ranging  between  -8192  and  + 8191.  Each 
data  record  contains  4096  words  or  8192  samples.  Since  each  data 
window  is  4096  samples  long,  each  record  represents  two  segments. 

The  computer  generated  processes  are  formed  from  the  FORTRAN 
routine  RAN  [431  which  produces  pseudo  random  numbers  distributed 
uniformly  on  the  interval  [0,11.  Twelve  consecutive  random  numbers 
are_^added  to  form  one  sample  of  an  approximately  Gaussian  process. 
The  parameters  of  the  process  are  adjusted  according  to 

V - V.^wi  (X,  - M)  (C.  4a) 

where 

0 - (1/2)  - n!  (VN)  , (C.  4b) 

V - { 12cr* / N) ,/z.  (C.4c) 

N is  the  number  of  additions,  is  a uniform  random  variable  on 
[0,11,  and  Y is  (asymptotically)  a Gaussian  random  variable  with 
mean  - //  and  variance  ■*  er1.  For  N = 12  and  u « 0,  V - <?  and  N » 


1/2. 


1 


I 
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For  this  work,  a normal  process  with  //  * 0 and  <r2  ■ 1000  is 
used.  The  value  1000  represents  a compromise  between  a distribution 
with  frequent  values  greater  than  8192  or  smaller  than  1 (which 
would  be  quantized  to  0).  Other  processes  are  then  formed  from  this 
basic  process  by  appropriate  manipulation  (such  as  digital  filtering 
with  an  appropriate  linear  system,  e.g.,  figures  2.3  and  3.2). 

The  nonstationary  process  discussed  in  section  4.6  is  produced 
by  scaling  each  4096  sample  section  by  a randomly  chosen  constant 
(this  is  equivalent  to  the  time-varying  linear  system  described  in 
section  4.2  being  an  amplifier).  The  random  gains  are  aiso 
generated  from  the  routine  RAN.  In  this  case,  the  uniform  random 
variables  are  scaled  and  exponentiated  to  have  the  desired 
characteristics.  Specifically,  if  X is  uniform  on  [0,1],  then 

V * exp  (X- log  (3 / A)  + log(A)]  (C.5) 

has  a hyperbolic  distribution  on  [A,B]  (see  section  B.3). 

The  real  data  come  from  a variety  of  sources.  The  singing  of 
figures  4.1  - 4.3  were  ciigitized  directly  from  phonograph  records. 
The  female  singer  "represented"  In- f i'guFe-’4.'4‘  ariS  “ The~sTrtPTg~ g'rtS’efltb'f e~ 
of  figure  4.5  were  digitized  directly  from  a live  microphone  and, 
hence,  are  free  from  the  noise  and  distortion  introduced  by 
intermediate  recording.  In  all  cases,  the  data  was  digitized  under 
the  direction  of  T.  G.  Stockham,  Jr. 

In  most  cases  the  dcta  was  sampled  at  10,000  Hertz  and  filtered 
at  4,000  Hertz.  The  sharp  cutoff  of  this  low-pass  filtering  is 

quite  evident  in  the  log  spectral  estimates.  In  the  case  of  the 
live  recordings,  sampling  was  at  37,500  Hertz  and  filtering  at 
15,000  Hertz.  Sampling  was  accomplished  with  a 14-bit  A/D  converter 
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interfaced  directly  with  the  computer. 

C.3  Hardware  and  Software  Description 

The  computations  in  this  document  were  performed  on  two  Digital 
Equipment  Corporation  PDP-10  computers;  one  a time-sharing  system 
with  a 262,144  word  memory  and  the  other  a single-user  system  with  a 
65,536  word  memory.  These  machines  have  floating  multiply  and 
floating  add  times  of  about  5 and  11  microseconds,  respectively 
[44]  . Mass  storage  o-:  data  was  on  high-speed  magnetic  disk; 
estimates  were  stored  on  magnetic  tape. 

The  log  spectral  estimates  were  computed  using  software  written 
in  FORTRAN  (the  author  made  extensive  modifications  of  original 
software  written  by  T.  G.  Stockham,  Jr.  for  implementation  of  the 
homomorphic  deconvolution  algorithm).  It  has  the  capability  of 
computing  both  the  average  log  and  iog  average  spectra  as  well  as 
displaying  intermediate  results.  It  also  allows  various  statistical 
parameters  to  be  computed  at  different  stages.  The  routine  makes 
use  of  several  subroutines  written  in  assembly  language  by  the 
programming  staff  of  the  Sensory  Information  Group  at  the  University 
of  Utah. 

Digital  filtering  was  accomplished  by  implementation  of  a 
high-speed  convolution  algorithm  [321.  This  approach  uses  the  fact 
that  convolution  is  maoped  into  multiplication  by  the  Fourier 
transform.  By  proper  bockkeeping  and  data  segmentation,  convolution 
is  realized  by  multiplying  the  Fourier  transforms  of  the  system 
impulse  response  and  each  data  segment,  inverse  transforming,  and 


summing. 


observed  for  filtering,  via  high-speed  convolution,  a process  of  the 
same  length  (42.0  minutes). 
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